Técnicas en Aprendizaje Estadístico - Introduction to Statistical Learning: Ejercicios Aplicados
#Librerías
library(tidyverse) # Manipulación, limpieza y gráficos
library(ISLR) # Bases de datos
library(plotly) # Interactividad gráfica
library(ggthemes) # Presentación
library(GGally) # Correlación
library(yardstick) # Métricas
library(MASS) # Clasificadores
library(class) # Clasificadores
library(tree) #CART - ISLR
library(broom) # Resultados tidy de modelos
library(latex2exp) # Expresiones LaTex en las gráficas
library(randomForest) # Libreria para bosques aleatorios
library(gbm) # Boosting
library(pls) # Principal component regression
library(ggforce) # Características adicionales para ggplot
library(e1071) # máquinas de soporte vectorial
library(caret) # Some models and plotsSección 4.7 - Ejercicio 10
a) Análisis descriptivo
El conjunto de datos Weekly contiene porcentajes de retorno semanales del índice S&P entre los años 1990 y 2010 con 1089 observaciones.
Year: es el año en el que fue obtenida la observación.Lag 1-5: Retornos porcentuales de las 5 semanas anteriores, respectivamente.Volume: Volumen de acciones intercambiadas (Número promedio de acciones diarias intercambiadas en billones).Today: Retorno porcentual de la semana.Direction: Indica si el mercado tuvo un retorno positivo o negativo en esta semana.
Encabezados y primeras 6 observaciones:
## # A tibble: 6 x 9
## Year Lag1 Lag2 Lag3 Lag4 Lag5 Volume Today Direction
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <fct>
## 1 1990 0.816 1.57 -3.94 -0.229 -3.48 0.155 -0.27 Down
## 2 1990 -0.27 0.816 1.57 -3.94 -0.229 0.149 -2.58 Down
## 3 1990 -2.58 -0.27 0.816 1.57 -3.94 0.160 3.51 Up
## 4 1990 3.51 -2.58 -0.27 0.816 1.57 0.162 0.712 Up
## 5 1990 0.712 3.51 -2.58 -0.27 0.816 0.154 1.18 Up
## 6 1990 1.18 0.712 3.51 -2.58 -0.27 0.154 -1.37 Down
Volumen promedio de acciones diarias intercambiadas en el tiempo
Se evidencia un crecimiento constante con un alza grande en el año 2005 y una estabilización de la trayectoria alrededor del año 2007. Se añadió un ruido a la posición de cada punto para poder evidenciar la densidad, además de transparencia. Se utilizó el método mgcv::gam para ajustar la tendencia.
p <- weekly %>%
ggplot(aes(x = Year, y = Volume)) +
geom_jitter(alpha = 0.6) +
geom_smooth() +
labs(x = NULL,
y = "Volumen (Billones)") +
theme_economist()
p_box <- weekly %>%
ggplot(aes(x = as_factor(Year), y = Volume, color = Today)) +
geom_boxplot() +
theme_economist() +
theme(axis.text.x = element_text(angle = 45)) +
labs(x = NULL,
y = "Volumen (Billones)")
ggplotly(p, width = 800, height = 500)El incremento en la volatilidad del índice acompañó su crecimiento a partir del 2007.
Correlación entre las variables
No se encontraron relaciones lineales claras entre las variables. Se presentan comportamientos esperados de la variable Today con la variable categórica Direction representada en el color, donde los valores negativos están asociados a un decrecimiento.
b) ¿Podemos predecir si el índice subirá o no mediante regresión logística?
Utilicemos regresión logística para predecir si habrá un retorno positivo del mercado en esta semana. Para ello, se utilizará el retorno porcentual de las 5 últimas semanas y el volumen sin separar un conjunto de entrenamiento y de prueba (Modelo sobreajustado).
market_logistic <- glm(Direction ~ Lag1 + Lag2 + Lag3 + Lag4 + Lag5 + Volume, data = weekly, family = binomial)
summary(market_logistic)##
## Call:
## glm(formula = Direction ~ Lag1 + Lag2 + Lag3 + Lag4 + Lag5 +
## Volume, family = binomial, data = weekly)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.6949 -1.2565 0.9913 1.0849 1.4579
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.26686 0.08593 3.106 0.0019 **
## Lag1 -0.04127 0.02641 -1.563 0.1181
## Lag2 0.05844 0.02686 2.175 0.0296 *
## Lag3 -0.01606 0.02666 -0.602 0.5469
## Lag4 -0.02779 0.02646 -1.050 0.2937
## Lag5 -0.01447 0.02638 -0.549 0.5833
## Volume -0.02274 0.03690 -0.616 0.5377
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1496.2 on 1088 degrees of freedom
## Residual deviance: 1486.4 on 1082 degrees of freedom
## AIC: 1500.4
##
## Number of Fisher Scoring iterations: 4
Si realizamos una prueba de hipótesis formal para la significancia de los parámetros, encontraremos que la probabilidad de rechazar la hipótesis nula es muy grande para todas las variables predictoras excepto para el Lag2 donde al parecer hay significancia. El valor en el tiempo \(t-2\) parece tener una relación con el tiempo \(t\).
c) Matriz de confusión de la regresión logística y calidad de la predicción del modelo sobreajustado
## Up
## Down 0
## Up 1
market_logistic_class <- tibble(value = predict(market_logistic, type = "response"),
pred = fct_relevel(as_factor(case_when(value > 0.5 ~ "Up",
TRUE ~ "Down")), "Down", "Up"),
real = weekly$Direction)
market_logistic_class %>%
conf_mat(truth = real, estimate = pred) %>%
autoplot(type = "heatmap") +
labs(title = "Matriz de confusión para el modelo de regresión logística")## # A tibble: 2 x 3
## .metric .estimator .estimate
## <chr> <chr> <dbl>
## 1 accuracy binary 0.561
## 2 kap binary 0.0350
Se obtiene una tasa de clasificación correcta accuracy del 56.1% utilizando regresión logística. Es importante interpretar este resultado correctamente, ya que la predicción obtenida indica que el mercado subió 557 días y bajó 54 días, y que no se está prediciendo el valor en el tiempo \(t+1\).
Cada predicción individual se debe interpretar de la siguiente manera: “Dadas las 5 semanas anteriores y el volumen de ésta, el mercado tendrá un retorno positivo”. La mayor cantidad de errores se cometió prediciendo que el mercado iba a subir y en realidad bajó, esta situación ocurrió 430 veces.
d) Modelo de regresión logística sólo con Lag2 y conjunto de prueba
weekly_train <- weekly %>%
filter(Year <= 2008)
weekly_test <- weekly %>%
filter(Year > 2008)
lag2_logistic <- glm(Direction ~ Lag2, data = weekly_train, family = "binomial")
summary(lag2_logistic)##
## Call:
## glm(formula = Direction ~ Lag2, family = "binomial", data = weekly_train)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.536 -1.264 1.021 1.091 1.368
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.20326 0.06428 3.162 0.00157 **
## Lag2 0.05810 0.02870 2.024 0.04298 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1354.7 on 984 degrees of freedom
## Residual deviance: 1350.5 on 983 degrees of freedom
## AIC: 1354.5
##
## Number of Fisher Scoring iterations: 4
lag2_class <- tibble(value = predict(lag2_logistic, type = "response", newdata = weekly_test),
pred = fct_relevel(as_factor(case_when(value > 0.5 ~ "Up",
TRUE ~ "Down")), "Down", "Up"),
real = weekly_test$Direction)
lag2_class %>%
conf_mat(truth = real, estimate = pred) %>%
autoplot(type = "heatmap") +
labs(title = "Matriz de confusión para la regresión logística con Lag2 como único predictor")## # A tibble: 2 x 3
## .metric .estimator .estimate
## <chr> <chr> <dbl>
## 1 accuracy binary 0.625
## 2 kap binary 0.141
Al igual que el modelo anterior, se tienen fallas al intentar predecir las bajas del mercado, sin embargo, este se validó en un conjunto de prueba de 104 observaciones desconocidas para el modelo. El accuracy en este caso es del 62.5%. Los resultados son sorpresivamente mejores que en el modelo sobreajustado del literal b).
d) Modelo LDA (Linear Discriminant Analysis) para predecir la tendencia del mercado
## Call:
## lda(Direction ~ Lag2, data = weekly_train)
##
## Prior probabilities of groups:
## Down Up
## 0.4477157 0.5522843
##
## Group means:
## Lag2
## Down -0.03568254
## Up 0.26036581
##
## Coefficients of linear discriminants:
## LD1
## Lag2 0.4414162
lag2_lda_class <- tibble(pred = predict(lag2_lda, newdata = weekly_test)$class,
real = weekly_test$Direction)
lag2_lda_class %>%
conf_mat(truth = real, estimate = pred) %>%
autoplot(type = "heatmap") +
labs(title = "Matriz de confusión para el modelo LDA con Lag2 como único predictor")## # A tibble: 2 x 3
## .metric .estimator .estimate
## <chr> <chr> <dbl>
## 1 accuracy binary 0.625
## 2 kap binary 0.141
Se obtienen las mismas clasificaciones para el modelo LDA y el modelo de regresión logística. Como el coeficiente discriminante lineal es Lag2 = 0.4414, es de esperar que en ambas categorías la gráfica de éstos sea similar. La matriz de confusión es igual a la del modelo logístico.
f) Modelo QDA (Quadratic Discriminant Analysis) para predecir la tendencia del mercado
## Call:
## qda(Direction ~ Lag2, data = weekly_train)
##
## Prior probabilities of groups:
## Down Up
## 0.4477157 0.5522843
##
## Group means:
## Lag2
## Down -0.03568254
## Up 0.26036581
lag2_qda_class <- tibble(pred = predict(lag2_qda, newdata = weekly_test)$class,
real = weekly_test$Direction)
lag2_qda_class %>%
conf_mat(truth = real, estimate = pred) %>%
autoplot(type = "heatmap") +
labs(title = "Matriz de confusión para el modelo QDA con Lag2 como único predictor")## # A tibble: 2 x 3
## .metric .estimator .estimate
## <chr> <chr> <dbl>
## 1 accuracy binary 0.587
## 2 kap binary 0
El clasificador QDA clasificó todos los valores como “Up”, si bien es una estrategia que otorga resultados buenos en términos porcentuales con respecto al accuracy, se espera que el modelo proponga en algunos casos predicciones de baja del mercado.
g) Modelo KNN (K-Nearest Neighbors) para predecir la tendencia del mercado
set.seed(1)
lag2_knn <- knn(train = as.matrix(weekly_train$Lag2), test = as.matrix(weekly_test$Lag2), cl = weekly_train$Direction, k = 3)
lag2_knn_class <- tibble(pred = lag2_knn,
real = weekly_test$Direction)
lag2_knn_class %>%
conf_mat(truth = real, estimate = pred) %>%
autoplot(type = "heatmap") +
labs(title = "Matriz de confusión para el modelo KNN con Lag2 como único predictor")## # A tibble: 2 x 3
## .metric .estimator .estimate
## <chr> <chr> <dbl>
## 1 accuracy binary 0.548
## 2 kap binary 0.0453
El clasificador \(KNN\) tiene una tasa de clasificación correcta del 50% utilizando \(k = 1\), sin embargo, al utilizar \(k = 3\) se obtiene un 54.8%. Este clasificador no tuvo inconvenientes en proponer predicciones variadas según el valor anterior del mercado.
h) ¿Qué métodos producen los mejores resultados para predecir la tendencia del mercado utilizando Lag2?
Los métodos de regresión logística y análisis de discriminante lineal obtuvieron los mejores resultados, con un accuracy del 62.5%. Se decide favorecer al modelo de regresión logística dada la facilidad de interpretación de los coeficientes y las facilidades gráficas de comunicación que permite la curva logística.
g <- augment_columns(lag2_logistic, type.predict = "response", newdata = weekly_test) %>%
arrange(.fitted) %>%
rownames_to_column() %>%
mutate(rowname = as.integer(rowname)) %>%
inner_join(lag2_class, by = c(.fitted = "value")) %>%
ggplot(aes(x = rowname, y= .fitted, color = Direction, label = Year, label1 = pred)) +
geom_point(alpha = 0.5, shape = 1, stroke = 2) +
geom_hline(yintercept = 0.5) +
theme_light() +
xlab("Índice") +
ylab("Probabilidad de que el mercado suba")
ggplotly(g, width = 800, height = 500)i) Mejora del mejor modelo y variables adicionales de interés
Selección Step-wise con la función stepAIC sin interacciones
full_logistic <- glm(Direction ~ .-Year-Today, data = weekly_train, family = "binomial")
stepAIC(full_logistic, direction = "both", trace = FALSE)##
## Call: glm(formula = Direction ~ Lag1 + Lag2, family = "binomial", data = weekly_train)
##
## Coefficients:
## (Intercept) Lag1 Lag2
## 0.21109 -0.05421 0.05384
##
## Degrees of Freedom: 984 Total (i.e. Null); 982 Residual
## Null Deviance: 1355
## Residual Deviance: 1347 AIC: 1353
step_logistic <- glm(Direction ~ Lag1 + Lag2, data = weekly_train, family = "binomial")
step_class <- tibble(value = predict(step_logistic, type = "response", newdata = weekly_test),
pred = fct_relevel(as_factor(case_when(value > 0.5 ~ "Up",
TRUE ~ "Down")), "Down", "Up"),
real = weekly_test$Direction)
step_class %>%
conf_mat(truth = real, estimate = pred) %>%
autoplot(type = "heatmap") +
labs(title = "Matriz de confusión para la regresión logística")## # A tibble: 2 x 3
## .metric .estimator .estimate
## <chr> <chr> <dbl>
## 1 accuracy binary 0.577
## 2 kap binary 0.0350
El modelo seleccionado por Stepwise selection tiene un accuracy menor que el de un solo predictor, sin embargo, su AIC es el mismo. El AIC evalúa únicamente el ajuste, por lo que se recomendaría el modelo de una sola variable predictora. Es importante resaltar que el modelo de menor AIC podría ser mejor en un conjunto de prueba diferente.
Mejor modelo por tanteo de interacciones
Se proponen variables con interacciones y se selecciona el modelo con mayor accuracy en el conjunto de prueba, en este caso utilizando la variable Lag2 y la interacción entre Lag1, Lag 4 y Volume
final_logistic <- glm(Direction ~ Lag2 + Lag1:Lag4:Volume, data = weekly_train, family = "binomial")
summary(final_logistic)##
## Call:
## glm(formula = Direction ~ Lag2 + Lag1:Lag4:Volume, family = "binomial",
## data = weekly_train)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.493 -1.265 1.023 1.089 1.446
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.2047966 0.0643540 3.182 0.00146 **
## Lag2 0.0553928 0.0291831 1.898 0.05768 .
## Lag1:Lag4:Volume 0.0007316 0.0014747 0.496 0.61982
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1354.7 on 984 degrees of freedom
## Residual deviance: 1350.3 on 982 degrees of freedom
## AIC: 1356.3
##
## Number of Fisher Scoring iterations: 4
final_class <- tibble(value = predict(final_logistic, type = "response", newdata = weekly_test),
pred = fct_relevel(as_factor(case_when(value > 0.5 ~ "Up",
TRUE ~ "Down")), "Down", "Up"),
real = weekly_test$Direction)
final_class %>%
conf_mat(truth = real, estimate = pred) %>%
autoplot(type = "heatmap") +
labs(title = "Matriz de confusión para la regresión logística")## # A tibble: 2 x 3
## .metric .estimator .estimate
## <chr> <chr> <dbl>
## 1 accuracy binary 0.644
## 2 kap binary 0.179
Con base en los resultados anteriores se selecciona el modelo por tanteo por su accuracy superior, sin embargo, esto puede ser un sobreajuste del conjunto de prueba; por ello, si este modelo fuera a producción se recomienda realizar validación cruzada para identificar si la interacción de variables añadida es significativa. Si no es lo es, el modelo seleccionado por Stepwise selection se considera apropiado.
Sección 4.7 - Ejercicio 11
En este apartado, desarrollaremos un modelo para predecir cuando un carro tiene alto o bajo rendimiento de combustible respecto al kilometraje.
Inicialmente observemos las características de este conjunto de datos. El conjunto de datos Auto contiene información sobre 392 vehículos cuyas variables son:
mpg: millas por galón.cylinders: número de cilindros (entre 4 y 8)displacement: desplazamiento del motor en pies.horsepower: caballos de fuerza del motor.weight: peso del vehículo en libras.acceleration: tiempo que tarda en acelerar de 0 a 60 medido en segundos.year: modelo del auto (módulo 100)origin: origen del vehículo.name: nombre del vehículo.
Se mostrarán los primeros 6 datos junto con su encabezado.
## # A tibble: 6 x 9
## mpg cylinders displacement horsepower weight acceleration year origin name
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <fct>
## 1 18 8 307 130 3504 12 70 1 chev~
## 2 15 8 350 165 3693 11.5 70 1 buic~
## 3 18 8 318 150 3436 11 70 1 plym~
## 4 16 8 304 150 3433 12 70 1 amc ~
## 5 17 8 302 140 3449 10.5 70 1 ford~
## 6 15 8 429 198 4341 10 70 1 ford~
a) Creación de una variable respuesta binaria
A continuación crearemos la variable mpg01 que tomará el valor de 1 si la variable mpg está por encima de la mediana, y 0 en caso contrario.
median_mpg <- median(auto$mpg)
auto <- auto %>% mutate(
mpg01 = case_when(
mpg >= median_mpg ~ 1,
mpg < median_mpg ~ 0
),
mpg01 = as.factor(mpg01),
origin = as.factor(origin)
)
# Retiramos la variable mpg
auto <- auto %>% dplyr::select(-mpg)
head(auto)## # A tibble: 6 x 9
## cylinders displacement horsepower weight acceleration year origin name mpg01
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <fct> <fct> <fct>
## 1 8 307 130 3504 12 70 1 chev~ 0
## 2 8 350 165 3693 11.5 70 1 buic~ 0
## 3 8 318 150 3436 11 70 1 plym~ 0
## 4 8 304 150 3433 12 70 1 amc ~ 0
## 5 8 302 140 3449 10.5 70 1 ford~ 0
## 6 8 429 198 4341 10 70 1 ford~ 0
b) Relaciones entre las variables
A continuación realicemos un análisis descriptivo para observar el comportamiento de las variables regresoras respecto a la variable mpg01
ggpairs(data = auto, columns = 1:6, mapping = aes(color = mpg01), legend = 1, lower = list(continuous = wrap("points", alpha = 0.6))) +
theme_light() +
theme(legend.position = "bottom")Por otro lado, demos una visualización para las variable origin, la cual es categórica.
ggpairs(data = auto, columns = c(7, 9), mapping = aes(color = mpg01), legend = 1, lower = list(continuous = wrap("points", alpha = 0.6))) +
theme_light() +
theme(legend.position = "bottom")Tomando en cuenta las gráficas, podemos afirmar que las variables cylinders, displacement, horsepower y weight son buenas variables cuantitativas para realizar una clasificación de la variable mpg01, pues al graficar discriminando por clases se puede observar que hay una partición notoria entre los carros de alto consumo y bajo consumo respecto a estas 4 variables. Note también que la variabele origen puede ser una buena variable predictora para mpg01 cuando el tipo de carro es americano.
c) División de conjuntos de entrenamiento y validación
A continuación dividiremos el conjunto Auto en dos conjuntos, uno de validación y otro de entrenamiento dejando el 75% de los datos para entrenar los modelos.
split_ratio <- 0.75
smp_size <- floor(split_ratio * nrow(auto))
train_ind <- sample(seq_len(nrow(auto)), size = smp_size)
auto_train <- auto[train_ind, ]
auto_test <- auto[-train_ind, ]
dim(auto_train)## [1] 294 9
## [1] 98 9
d) Modelo LDA
A continuación realizaremos un modelo LDA considerando las varables que fueron extraídas en el análisis descriptivo y que se consideraron importantes a la hora de realizar una clasificación para la variable respuesta mpg01.
lda_model = lda(mpg01 ~ cylinders + displacement + horsepower + weight + origin, data = auto_train)
lda_model## Call:
## lda(mpg01 ~ cylinders + displacement + horsepower + weight +
## origin, data = auto_train)
##
## Prior probabilities of groups:
## 0 1
## 0.5204082 0.4795918
##
## Group means:
## cylinders displacement horsepower weight origin2 origin3
## 0 6.745098 270.0784 128.03922 3595.412 0.05882353 0.04575163
## 1 4.163121 115.5887 79.26241 2339.532 0.26241135 0.38297872
##
## Coefficients of linear discriminants:
## LD1
## cylinders -0.4317395853
## displacement -0.0018236085
## horsepower 0.0018307395
## weight -0.0007360167
## origin2 0.3722004164
## origin3 0.5312690065
Grafiquemos a continuación el resultado del modelo.
Ahora construyamos la matriz de confusión para el modelo con el conjunto de datos de validación para observar el desempeño.
mpg_lda_class <- tibble(pred = predict(lda_model, auto_test %>% dplyr::select(-mpg01))$class ,
real = auto_test$mpg01)
mpg_lda_class %>% conf_mat(truth = real, estimate = pred) %>%
autoplot(type = "heatmap") +
labs(title = "Matriz de confusión para el modelo LDA")## # A tibble: 2 x 3
## .metric .estimator .estimate
## <chr> <chr> <dbl>
## 1 accuracy binary 0.867
## 2 kap binary 0.729
Note entonces que la tasa de clasificación correcta es de un 86.73%, lo que muestra que el modelo tiene una buena calidad.
e) Modelo QDA
Con el fin de realizar una comparativa con el modelo anterior, realicemos un modelo QDA para predecir la varibale mpg01 usando las variables más significativas encontradas en el análisis descriptivo.
qda_model = qda(mpg01 ~ cylinders + displacement + horsepower + weight + origin, data = auto_train)
qda_model## Call:
## qda(mpg01 ~ cylinders + displacement + horsepower + weight +
## origin, data = auto_train)
##
## Prior probabilities of groups:
## 0 1
## 0.5204082 0.4795918
##
## Group means:
## cylinders displacement horsepower weight origin2 origin3
## 0 6.745098 270.0784 128.03922 3595.412 0.05882353 0.04575163
## 1 4.163121 115.5887 79.26241 2339.532 0.26241135 0.38297872
Ahora construyamos la matriz de confusión para el modelo QDA con el conjunto de datos de validación.
mpg_qda_class <- tibble(pred = predict(qda_model, auto_test %>% dplyr::select(-mpg01))$class ,
real = auto_test$mpg01)
mpg_qda_class %>% conf_mat(truth = real, estimate = pred) %>%
autoplot(type = "heatmap") +
labs(title = "Matriz de confusión para el modelo QDA")## # A tibble: 2 x 3
## .metric .estimator .estimate
## <chr> <chr> <dbl>
## 1 accuracy binary 0.867
## 2 kap binary 0.729
Note entonces que la tasa de clasificación correcta es de un 86.73%, y en comparación con el modelo LDA. Se puede decir que tienen un desempeño muy parecido de acuerdo a la tasa de clasificación correcta.
f) Modelo de regresión logística
Por otra parte, construiremos un modelo de regresión logística para realizar una predicción de la variable mpg01.
logistic_model <- glm(mpg01 ~ cylinders + displacement + horsepower + weight + origin, data = auto_train, family = binomial)
summary(logistic_model)##
## Call:
## glm(formula = mpg01 ~ cylinders + displacement + horsepower +
## weight + origin, family = binomial, data = auto_train)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.53196 -0.25902 -0.00542 0.33608 3.03362
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 11.9583896 2.1576414 5.542 2.98e-08 ***
## cylinders -0.1003094 0.4308485 -0.233 0.81590
## displacement -0.0117784 0.0122590 -0.961 0.33665
## horsepower -0.0546728 0.0176520 -3.097 0.00195 **
## weight -0.0015796 0.0009063 -1.743 0.08135 .
## origin2 0.1646890 0.6822990 0.241 0.80927
## origin3 0.6067756 0.6915181 0.877 0.38024
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 407.08 on 293 degrees of freedom
## Residual deviance: 152.08 on 287 degrees of freedom
## AIC: 166.08
##
## Number of Fisher Scoring iterations: 7
Calculemos de igual manera la matriz de confusión para realizar una comparación con los métodos previos usando los datos de entrenamiento.
predictions_prob_glm <- predict(logistic_model, auto_test %>% dplyr::select(-mpg01), type = "response")
predictions_glm <- rep(0, nrow(auto_test))
predictions_glm[predictions_prob_glm > 0.5] <- 1
mpg_logistic_class <- tibble(pred = as.factor(predictions_glm),
real = auto_test$mpg01)
mpg_logistic_class %>% conf_mat(truth = real, estimate = pred) %>%
autoplot(type = "heatmap") +
labs(title = "Matriz de confusión para el modelo de regresión logística")test_results_logistic <- mpg_logistic_class %>%
metrics(truth = real, estimate = pred)
test_results_logistic## # A tibble: 2 x 3
## .metric .estimator .estimate
## <chr> <chr> <dbl>
## 1 accuracy binary 0.888
## 2 kap binary 0.772
De acuerdo a este resultado, se puede ver que todos los métodos tienen un puntaje muy similar entre ellos pues en este caso, se puede observar una tasa de clasificación correcta del 88.78%.
g) Modelo \(KNN\)
Realicemos a continuación un modelo de \(KNN\) para predecir la variable mpg, en este caso, debemos probar con varios casos de \(K\) para hallar un parámetro óptimo. Para esto, realizaremos una función que reciba como parámetro un numero \(k\) y esta función retornará la tasa de clasificación correcta para el \(K\) dado, realizando la clasificación sobre el conjunto de prueba.
get_accuracy_knn <- function(k) {
set.seed(10)
knn_model <- knn(auto_train %>% dplyr::select(-mpg01, -name),
auto_test %>% dplyr::select(-mpg01, -name),
auto_train$mpg01,
k)
knn_class <- tibble(pred = as.factor(knn_model), real = auto_test$mpg01)
results_knn <- knn_class %>%
metrics(truth = real, estimate = pred)
return(results_knn$.estimate[1])
}Ahora tomaremos varios valores de \(K\) en un rango entre 1 y 60, y para estos valores, computaremos su tasa de clasificación correcta y los graficaremos para observar cual es el mejor valor de \(k\) en este rango.
k_values = 1:60
accuracy_results <- tibble(k = k_values, accuracy = sapply(k_values, get_accuracy_knn))
g <- ggplot(data = accuracy_results) +
geom_line(mapping = aes(x = k, y = accuracy)) +
xlab("K") +
ylab("Tasa de clasificación correcta") +
theme_light()
ggplotly(g, width = 800, height = 500)Una vez observada la gráfica, encontraremos el mejor valor de \(k\) que mejora la tasa de clasificación correcta, y además para este valor de \(k\) realizaremos la matriz de confusión para ver el desempeño del modelo de manera global.
set.seed(10)
best_k <- min(
accuracy_results[which(accuracy_results$accuracy == max(accuracy_results$accuracy)), "k"]
)
knn_best_model <- knn(auto_train %>% dplyr::select(-mpg01, -name),
auto_test %>% dplyr::select(-mpg01, -name),
auto_train$mpg01,
best_k)
mpg_knn_class <- tibble(pred = as.factor(knn_best_model), real = auto_test$mpg01)
test_results_knn <- mpg_knn_class %>%
metrics(truth = real, estimate = pred)
mpg_knn_class %>% conf_mat(truth = real, estimate = pred) %>%
autoplot(type = "heatmap") +
labs(title = "Matriz de confusión para el modelo KNN")Concluimos entonces que para los valores de \(k\) probados, la mejor tasa de clasificación correcta fue de 89.8%, la cuál sigue siendo muy similar a los puntajes anteriores y desde un punto de vista global, estos puntajes denotan que los modelos tienen una muy buena calidad.
Sección 4.7 - Ejercicio 12
a) Implementación de la función Power()
A continuación implementaremos una función Power() que calculará la potencia \(2^3\)
b) Implementación de la función Power2() generalizada
Ahora, implementaremos la función Power2() que recibe dos parámetros x y a y computará \(x^a\)
c) Usos de la funcion Power2()
En el siguiente bloque de código calcularemos \(10^3\).
## [1] 1000
Ahora calcularemos el valor de \(8^{17}\).
## [1] 2.2518e+15
Y finalmente calcularemos el valor de \(131^3\).
## [1] 2248091
d) Implementación de la función Power3()
Ahora implementaremos una versión mucho más general para realizar las potencias que esta basada en la función Power2(), pero en este caso la función no imprimirá el resultado sino que retornara dicho valor como un objeto de R
e) Usando Plot3() para graficar la función \(f(x) = x ^ 2\)
En este bloque de código, graficaremos \(f(x) = x ^ 2\) usando una escala logarítmica para el eje \(y\)
ggplot(data = tibble(x = 1:10, y = sapply(1:10, Power3, 2))) +
geom_line(mapping = aes(x = x, y = y)) +
xlab("x") +
ylab(TeX("x^2")) +
scale_y_log10()f) Implementación de una forma generalizada para graficación
A continuación implementaremos un código para graficar la función \(f(x) = x ^ a\) tal que \(x \in I\) para un intervalo \(I \subseteq \mathbb{R}\).
PlotPower <- function(interval, a) {
ggplot(data = tibble(x = interval, y = sapply(interval, Power3, a))) +
geom_line(mapping = aes(x = x, y = y)) +
xlab("x") +
ylab(TeX(paste0("x^", toString(a)))) +
labs(title = "Gráfica de la función y = f(c)")
}
PlotPower(-10:10, 7)Sección 4.7 - Ejercicio 13
Utilizaremos el conjunto de datos Boston para predecir si un suburbio tiene un porcentaje de crimen mayor o menor que la media. El conjunto de datos tiene las siguientes variables:
crim: Ratio de crimen per cápita por pueblo.zn: Proporción de zonas residenciales asignadas para lotes mayores a 25.000 sq.ft.indus: Proporción de negocios no comerciales en acres por pueblo.chas: 1 si está en trayecto del río Charles y 0 en otro caso.nox: Concentración de óxidos de nitrógeno (partes por millón).rm: Número promedio de cuartos por hogar.age: Proporción de unidades ocupadas construidas antes de 1940.dis: Media ponderada de las distancias a 5 de los centros de empleo de Boston.rad: Índice de accesibilidad a autopistas periféricas.tax: Valor total de los impuestos por cada $10.000ptratio: Ratio de pupilos-profesores por pueblo.black: \(1000(Bk - 0.63)^2\) dónde \(Bk\) es la proporción de negros por pueblo.lstat: Menor estatus de la población (porcentaje).medv: Valor medio de una propiedad ocupada en $1000s.
Creación de la variable de interés
boston <- as_tibble(Boston) %>%
mutate(crimen = factor(case_when(crim > median(crim) ~ "alto",
TRUE ~ "bajo")))La mediana de la variable crim es 0.25651, se creará la variable crimen con niveles “alto” si es mayor a la mediana y “bajo” en otro caso.
Análisis descriptivo
Las variables de edad e indus presentan comportamientos discriminantes de interés, en la variable edad, los valores mayores a 70 tienen una gran prevalencia en la categoría alto y menores a este valor en la categoría bajo, por lo que al incluirla en uno de los modelos de clasificación se espera que sea significativa.
Este mismo comportamiento se identifica en la variable indus, donde los valores menores a 15 parecen tener una gran prevalencia en la categoría bajo, y entre 18 y 22 se encuentra la mayor cantidad de alto.
crim_hist <- boston %>%
ggplot(aes(x = age, fill = crimen)) +
geom_histogram(color = "black") +
labs(title = "Proporción de unidades ocupadas construidas antes de 1940.")
ggplotly(crim_hist)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
indus_hist <- boston %>%
ggplot(aes(x = indus, fill = crimen)) +
geom_histogram(color = "black") +
labs(title = "Proporción de negocios no comerciales en acres por pueblo.")
ggplotly(indus_hist)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Regresión logística
Construyamos inicialmente un clasificador únicamente con las dos características identificadas en el análisis descriptivo para evaluar las suposiciones. Realicemos también la separación en los conjutos de prueba y de entrenamiento.
set.seed(1995)
boston_train <- boston %>% sample_frac(size = 0.8)
boston_test <- boston %>% anti_join(boston_train)## Joining, by = c("crim", "zn", "indus", "chas", "nox", "rm", "age", "dis", "rad", "tax", "ptratio", "black", "lstat", "medv", "crimen")
logistic_boston_manual <- glm(crimen ~ age + indus, data = boston_train, family = "binomial")
summary(logistic_boston_manual)##
## Call:
## glm(formula = crimen ~ age + indus, family = "binomial", data = boston_train)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.6774 -0.5407 -0.3414 0.5724 2.6981
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.727028 0.491097 9.625 < 2e-16 ***
## age -0.047380 0.007062 -6.709 1.96e-11 ***
## indus -0.131905 0.025355 -5.202 1.97e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 561.25 on 404 degrees of freedom
## Residual deviance: 337.58 on 402 degrees of freedom
## AIC: 343.58
##
## Number of Fisher Scoring iterations: 5
Efectivamente ambos parámetros son altamente significativos. evaluemos la calidad de las predicciones.
boston_manual_class <- tibble(value = predict(logistic_boston_manual, type = "response", newdata = boston_test),
pred = fct_relevel(as_factor(case_when(value > 0.5 ~ "bajo",
TRUE ~ "alto")), "alto", "bajo"),
real = boston_test$crimen)
boston_manual_class %>%
conf_mat(truth = real, estimate = pred) %>%
autoplot(type = "heatmap") +
labs(title = "Matriz de confusión para la regresión logística con age y indus")## # A tibble: 2 x 3
## .metric .estimator .estimate
## <chr> <chr> <dbl>
## 1 accuracy binary 0.772
## 2 kap binary 0.547
Los resultados son aceptables, sin embargo, evaluemos el modelo step-wise propuesto para identificar otras variables de interés.
boston_full_model <- glm(crimen ~ . -crimen - crim, data = boston_train, family = "binomial")
summary(boston_full_model)##
## Call:
## glm(formula = crimen ~ . - crimen - crim, family = "binomial",
## data = boston_train)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.4328 -0.0032 0.0000 0.1377 1.6830
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 27.787319 8.586662 3.236 0.00121 **
## zn 0.062198 0.033850 1.837 0.06614 .
## indus 0.068162 0.055431 1.230 0.21882
## chas -1.101626 0.829059 -1.329 0.18393
## nox -46.572894 8.258444 -5.639 1.71e-08 ***
## rm 0.079837 0.774442 0.103 0.91789
## age -0.033585 0.014415 -2.330 0.01981 *
## dis -0.659126 0.244784 -2.693 0.00709 **
## rad -0.663217 0.167363 -3.963 7.41e-05 ***
## tax 0.006860 0.003033 2.262 0.02372 *
## ptratio -0.458391 0.151742 -3.021 0.00252 **
## black 0.035793 0.014010 2.555 0.01062 *
## lstat -0.006730 0.058173 -0.116 0.90790
## medv -0.141170 0.073683 -1.916 0.05538 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 561.25 on 404 degrees of freedom
## Residual deviance: 160.46 on 391 degrees of freedom
## AIC: 188.46
##
## Number of Fisher Scoring iterations: 9
##
## Call: glm(formula = crimen ~ zn + nox + age + dis + rad + tax + ptratio +
## black + medv, family = "binomial", data = boston_train)
##
## Coefficients:
## (Intercept) zn nox age dis rad
## 26.279813 0.070380 -41.630025 -0.033514 -0.621363 -0.762505
## tax ptratio black medv
## 0.008532 -0.402343 0.031772 -0.134618
##
## Degrees of Freedom: 404 Total (i.e. Null); 395 Residual
## Null Deviance: 561.2
## Residual Deviance: 163.3 AIC: 183.3
boston_stepwise <- glm(formula = crimen ~ zn + nox + age + dis + rad + tax + ptratio + black + medv,
family = "binomial",
data = boston_train)
boston_stepwise_class <- tibble(value = predict(boston_stepwise, type = "response", newdata = boston_test),
pred = fct_relevel(as_factor(case_when(value > 0.5 ~ "bajo",
TRUE ~ "alto")), "alto", "bajo"),
real = boston_test$crimen)
boston_stepwise_class %>%
conf_mat(truth = real, estimate = pred) %>%
autoplot(type = "heatmap") +
labs(title = "Matriz de confusión para la regresión logística con el modelo stepwise")## # A tibble: 2 x 3
## .metric .estimator .estimate
## <chr> <chr> <dbl>
## 1 accuracy binary 0.901
## 2 kap binary 0.800
boston_plot <- augment_columns(boston_stepwise, type.predict = "response", newdata = boston_test) %>%
arrange(.fitted) %>%
rownames_to_column() %>%
mutate(rowname = as.integer(rowname)) %>%
inner_join(boston_stepwise_class, by = c(.fitted = "value")) %>%
ggplot(aes(x = rowname, y= .fitted, color = crimen, label1 = pred)) +
geom_point(alpha = 0.5, shape = 1, stroke = 2) +
geom_hline(yintercept = 0.5) +
theme_light() +
xlab("Índice") +
ylab("")
ggplotly(boston_plot, width = 800, height = 500)El modelo stepwise tiene una mejora considerable en las métricas al incluir otras variables adicionales, sin embargo, se destaca que el poder discrimatorio de las variables adicionales con respecto al modelo inicial propuesto sólo mejoró la predicción en un 12.9%. En un ambiente de producción se recomendaría evaluar el costo de utilizar estas variables en el modelo versus el beneficio de predecir correctamente.
Modelo LDA
boston_lda <- lda(crimen ~ zn + nox + age + dis + rad + tax + ptratio + black + medv, data = boston_train)
boston_lda## Call:
## lda(crimen ~ zn + nox + age + dis + rad + tax + ptratio + black +
## medv, data = boston_train)
##
## Prior probabilities of groups:
## alto bajo
## 0.5111111 0.4888889
##
## Group means:
## zn nox age dis rad tax ptratio black
## alto 1.178744 0.6392367 86.25942 2.528803 14.787440 507.3575 18.97874 325.7682
## bajo 22.209596 0.4703273 51.19040 5.149906 4.136364 306.2677 17.84192 390.5477
## medv
## alto 20.33720
## bajo 25.31162
##
## Coefficients of linear discriminants:
## LD1
## zn 0.005331755
## nox -8.025497066
## age -0.016986785
## dis -0.053038787
## rad -0.074033448
## tax 0.001113865
## ptratio -0.062819427
## black 0.001227959
## medv -0.033685626
boston_lda_class <- tibble(pred = predict(boston_lda, newdata = boston_test)$class,
real = boston_test$crimen)
boston_lda_class %>%
conf_mat(truth = real, estimate = pred) %>%
autoplot(type = "heatmap") +
labs(title = "Matriz de confusión para el modelo LDA con age e indus")## # A tibble: 2 x 3
## .metric .estimator .estimate
## <chr> <chr> <dbl>
## 1 accuracy binary 0.891
## 2 kap binary 0.778
Los resultados son muy similares a los del modelo de regresión logística siendo ligeramente peores, sin embargo, la gráfica construida por el LDA ofrece una interpretabilidad mayor a la obtenida en el literal 10.
Modelo KNN (K-Nearest Neighbors)
boston_train_scaled <- boston_train %>%
dplyr::select(-crim, -crimen, -indus, -chas, -lstat) %>%
scale()
boston_test_scaled <- boston_test %>%
dplyr::select(-crim, -crimen, -indus, -chas, -lstat) %>%
scale()
boston_knn <- knn(train = boston_train_scaled, test = boston_test_scaled, cl = boston_train$crimen, k = 5)
boston_knn_class <- tibble(pred = boston_knn,
real = boston_test$crimen)
boston_knn_class %>%
conf_mat(truth = real, estimate = pred) %>%
autoplot(type = "heatmap") +
labs(title = "Matriz de confusión para el modelo KNN con todos las variables")## # A tibble: 2 x 3
## .metric .estimator .estimate
## <chr> <chr> <dbl>
## 1 accuracy binary 0.931
## 2 kap binary 0.861
El clasificador \(KNN\) tiene una tasa de clasificación correcta del 91.1% utilizando \(k = 3\) y todas las variables, sin embargo, al utilizar \(k = 5\) y las mismas variables que la regresión logística se obtiene un 93.1% de predicción correcta.
Sección 8.4 - Ejercicio 7
En este literal, tomaremos el conjunto de datos Boston y analizaremos de una manera mas detallada el efecto de los parámetros en el \(MSE\).
Inicialmente creamos el conjunto de datos a analizar
## # A tibble: 6 x 14
## crim zn indus chas nox rm age dis rad tax ptratio black
## <dbl> <dbl> <dbl> <fct> <dbl> <dbl> <dbl> <dbl> <int> <dbl> <dbl> <dbl>
## 1 0.00632 18 2.31 0 0.538 6.58 65.2 4.09 1 296 15.3 397.
## 2 0.0273 0 7.07 0 0.469 6.42 78.9 4.97 2 242 17.8 397.
## 3 0.0273 0 7.07 0 0.469 7.18 61.1 4.97 2 242 17.8 393.
## 4 0.0324 0 2.18 0 0.458 7.00 45.8 6.06 3 222 18.7 395.
## 5 0.0690 0 2.18 0 0.458 7.15 54.2 6.06 3 222 18.7 397.
## 6 0.0298 0 2.18 0 0.458 6.43 58.7 6.06 3 222 18.7 394.
## # ... with 2 more variables: lstat <dbl>, medv <dbl>
Ahora graficaremos el error de prueba con una variedad de valores para mtry y ntree. Para calcular el error de prueba en cada uno de los valores, crearemos una función que evalua el error de prueba en cada caso.
mtry_values = 2:(ncol(boston) - 1)
ntree_values = 3:100
test_error_random_forest <- function(row, train, test) {
mtry_test <- row[1]
ntree_test <- row[2]
test_tree <- randomForest(formula = medv ~ .,
data = train,
mtry = mtry_test,
ntree = ntree_test)
test_predicted_medv <- predict(test_tree, newdata = test)
error <- mean((test_predicted_medv - test$medv) ^ 2)
return(error)
}Ahora dividiremos el conjunto de datos en datos de entrenamiento y datos de prueba. Se obtienen las siguientes dimensiones de los conjuntos de entrenamiento y prueba respectivamente.
split_ratio <- 0.75
smp_size <- floor(split_ratio * nrow(boston))
train_ind <- sample(seq_len(nrow(boston)), size = smp_size)
boston_train <- boston[train_ind, ]
boston_test <- boston[-train_ind, ]
dim(boston_train)## [1] 379 14
## [1] 127 14
Creamos un tabla para los valores a probar.
test_values <- expand.grid(mtry_values, ntree_values)
colnames(test_values) <- c("mtry_values", "ntree_values")
test_values <- as_tibble(test_values)
test_values["error"] <- apply(test_values, 1, test_error_random_forest, boston_train, boston_test)
head(test_values)## # A tibble: 6 x 3
## mtry_values ntree_values error
## <int> <int> <dbl>
## 1 2 3 26.5
## 2 3 3 23.5
## 3 4 3 16.9
## 4 5 3 21.2
## 5 6 3 16.9
## 6 7 3 15.2
Y finalmente graficamos los errores mediante una gráfica de contornos como se muestra a continuación.
fig <- plot_ly(data = test_values, x=~mtry_values,y=~ntree_values, z=~error, type = "contour", colorscale='Jet') %>%
layout(
xaxis = list(title = "Valores de mtry"),
yaxis = list(title = "Valores de ntree"))
fig %>% colorbar(title = "MSE")Observemos que para esta selección de conjuntos de entrenamiento y de validación, se cumple que el bagging no es el metodo mas efectivo, pues allí la grafica se muestra mas roja, y por tanto evidencia un valor mas alto del MSE. Para este caso, el valor de mtree más adecuado se encuentra entre 4 y 5. Note que los valores de ntree no tienen una relación significativa con el \(MSE\), pues los mejores valores del error de validación se pueden encontrar con valores de ntreecercanos a 100 pero también con valores cercanos a 50. Sin embargo, valores bajos del ntree dan errores muy altos sin importar el valor de mtree.
Sección 8.4 - Ejercicio 8
Analizaremos el conjunto de datos Carseats, el cuál es un conjunto de datos que posee información respecto a sillas para niños, las cuales que se ubican dentro de los automóviles, y realizaremos una regresión para la variable Sales. Las variables que tiene este conjunto de datos son las siguientes:
Sales: unidades vendidas en miles en cada una de las ubicaciones.CompPrice: precio por competidor en cada ubicación.Income: ingresos de la comunidad en miles de dólares.Advertising: presupuesto para propaganda local en cada ubicación.Population: población en cada región.Price: precio de los asientos en cada sitio.ShelveLoc: calidad de los puntos de venta en cada ubicación.Age: edad promedio de la población en cada ubicación.Education: nivel de educación en cada ubicación.Urban: determina si la tienda se encuentra en territorio rural o urbano.US: determina si la tienda esta en EE.UU.
## # A tibble: 6 x 11
## Sales CompPrice Income Advertising Population Price ShelveLoc Age Education
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <fct> <dbl> <dbl>
## 1 9.5 138 73 11 276 120 Bad 42 17
## 2 11.2 111 48 16 260 83 Good 65 10
## 3 10.1 113 35 10 269 80 Medium 59 12
## 4 7.4 117 100 4 466 97 Medium 55 14
## 5 4.15 141 64 3 340 128 Bad 38 13
## 6 10.8 124 113 13 501 72 Bad 78 16
## # ... with 2 more variables: Urban <fct>, US <fct>
a) División del conjunto de datos en entrenamiento y validación
A continuación dividiremos el conjunto de datos en un conjunto de entrenamiento y otro de prueba, obteniendo las siguientes dimensiones respectivamente.
split_ratio <- 0.75
smp_size <- floor(split_ratio * nrow(carseats))
train_ind <- sample(seq_len(nrow(carseats)), size = smp_size)
carseats_train <- carseats[train_ind, ]
carseats_test <- carseats[-train_ind, ]
dim(carseats_train)## [1] 300 11
## [1] 100 11
b) Entrenamiento de un árbol de regresión
Ahora, entrenaremos un árbol de regresión.
## node), split, n, deviance, yval
## * denotes terminal node
##
## 1) root 300 2383.00 7.412
## 2) ShelveLoc: Bad,Medium 238 1373.00 6.706
## 4) Price < 92.5 37 166.30 8.921
## 8) Income < 47.5 7 37.60 6.526 *
## 9) Income > 47.5 30 79.18 9.480 *
## 5) Price > 92.5 201 992.00 6.298
## 10) ShelveLoc: Bad 61 227.60 4.782
## 20) Price < 129.5 43 145.30 5.262
## 40) CompPrice < 125.5 26 64.85 4.477 *
## 41) CompPrice > 125.5 17 39.86 6.462 *
## 21) Price > 129.5 18 48.78 3.636 *
## 11) ShelveLoc: Medium 140 562.90 6.959
## 22) Advertising < 14.5 118 406.30 6.626
## 44) Price < 127 77 195.70 7.198
## 88) CompPrice < 124.5 42 70.21 6.453 *
## 89) CompPrice > 124.5 35 74.12 8.093 *
## 45) Price > 127 41 138.20 5.553
## 90) Age < 64.5 29 57.74 6.153 *
## 91) Age > 64.5 12 44.83 4.103 *
## 23) Advertising > 14.5 22 73.45 8.745
## 46) Age < 54.5 14 25.93 9.592 *
## 47) Age > 54.5 8 19.86 7.261 *
## 3) ShelveLoc: Good 62 436.80 10.120
## 6) Price < 109.5 19 59.60 12.420
## 12) CompPrice < 120.5 13 17.09 11.570 *
## 13) CompPrice > 120.5 6 12.62 14.270 *
## 7) Price > 109.5 43 231.80 9.102
## 14) Advertising < 13.5 35 150.70 8.500
## 28) Age < 60.5 21 70.78 9.500 *
## 29) Age > 60.5 14 27.42 7.000 *
## 15) Advertising > 13.5 8 12.84 11.740 *
La gráfica del arbol se muestra a continuación.
Con este arbol entrenado, calculemos el error de prueba como se muestra a continuación
test_predict <- predict(tree_carseats, carseats_test)
mse_test <- mean((test_predict - carseats_test$Sales) ^ 2)
mse_test## [1] 4.477323
Y también mostremos el \(RMSE\) para este modelo tomando el conjunto de entrenamiento.
## [1] 2.115969
Como conclusión se puede observar que una variable que determina el número de ventas de manera diferenciadora es la variable ShelveLoc, note entonces que cuando un local tiene buena calidad en las estanterías del punto de venta, la cantidad de ventas es mucho mayor respecto a locales con una calidad más baja. Aún así, los productos con un precio más bajo tienden a tener una cantidad de unidades vendidas mucho más alta. Por tanto, dos aspectos que favorecen a las ventas es tener un local con buena calidad y precios por debajo de 97.5 dólares. Finalemente, considerando que la raíz del error cuadrático medio es aproximadamente de 2.1159687, quiere decir que el modelo se equivoca en un radio de 2.1159687, lo que realmente da una predicción aceptable sobre los asientos vendidos.
c) Uso de validación cruzada para determinar el nivel de complejidad óptimo
En esta sección podaremos el árbol entrenado usando validación cruzada, con el objetivo de mejorar la calidad de la recomendación.
## $size
## [1] 16 15 14 13 11 10 9 8 7 6 5 4 3 2 1
##
## $dev
## [1] 1498.228 1531.632 1529.217 1513.865 1516.933 1500.161 1484.006 1491.751
## [9] 1547.555 1557.176 1556.349 1648.648 1873.875 1864.896 2394.920
##
## $k
## [1] -Inf 27.65922 29.89327 35.65000 37.04248 49.53619 51.34393
## [8] 52.48000 68.26306 72.40865 83.19008 145.38605 201.45179 214.92784
## [15] 573.23081
##
## $method
## [1] "deviance"
##
## attr(,"class")
## [1] "prune" "tree.sequence"
Grafiquemos ahora la sumatoria de errores cuadráticos respecto al tamaño del árbol.
g <- ggplot(data = data.frame(size = cv_prune_tree_carseats$size, dev = cv_prune_tree_carseats$dev)) +
geom_line(mapping = aes(x = size, y = dev)) +
geom_point(mapping = aes(x = size, y = dev)) +
labs(x = "Numero de nodos hoja", y = "Sumatoria de errores cuadráticos", title = "Error de predicción respecto al tamaño del árbol")
ggplotly(g)Note entonces que la validación cruzada arroja que el mejor árbol es el árbol original sin podar. Como se puede apreciar en la gráfica, el arbol mas grande es aquel que tiene mejor sumatoria de errores cuadráticos, por tanto no se recomienda podar el árbol.
d) Bagging e importancia de las variables
Ahora, usaremos bagging con el objetivo de observar si hay una mejora en la predicción y también se requiere observar cuáles son las variables mas importantes en este modelo.
Incialmente realizaremos bagging y compararemos el \(RMSE\) con los métodos previos.
bagging_carseats <- randomForest(Sales ~ ., data = carseats_train, importance = TRUE, mtry = ncol(carseats) - 1)
bagging_carseats##
## Call:
## randomForest(formula = Sales ~ ., data = carseats_train, importance = TRUE, mtry = ncol(carseats) - 1)
## Type of random forest: regression
## Number of trees: 500
## No. of variables tried at each split: 10
##
## Mean of squared residuals: 2.551421
## % Var explained: 67.88
Calculemos ahora el error de validación para este método
test_bagging_carseats <- predict(bagging_carseats, carseats_test)
test_mse_bagging_carseats <- mean((test_bagging_carseats - carseats_test$Sales) ^ 2)
test_mse_bagging_carseats## [1] 2.380399
Note entonces que el error de validación es mayor en este caso respecto al árbol de regresión. Ahora usemos el método importance() para determinar la importancia de las variables en este conjunto de datos.
## %IncMSE IncNodePurity
## CompPrice 30.4360564 223.022071
## Income 11.1409987 147.056363
## Advertising 21.6230477 183.453069
## Population -1.4934508 78.211847
## Price 62.9001945 676.637190
## ShelveLoc 71.8226799 754.037937
## Age 20.1073551 174.146019
## Education 5.5360020 71.905929
## Urban -0.4547497 9.617858
## US 4.7873814 10.535014
Veamos una grafica de la importancia para cada una de las variables.
Note entonces que el incremento en la impureza y el porcentaje de incremento del MSE revelan que las variables ShelveLoc y Price son las que tienen mas relevancia dentro del modelo, y esta relevancia es bastante marcada respecto a las demás variables.
e) Uso de bósques aleatorios
En esta ocasión usaremos bosques aleatorios para determinar si este método es más efectivo y además, observaremos cuales son las variables de mayor importancia para este método considerando el conjunto de datos Carseats.
Primero entrenemos un modelo con los parámetros por defecto, obteniendo el siguiente \(MSE\)
random_forest_carseats <- randomForest(Sales ~ ., data = carseats_train)
result_test_carseats <- predict(random_forest_carseats, newdata = carseats_test)
mean((result_test_carseats - carseats_test$Sales) ^ 2)## [1] 2.882621
Luego, analicemos la importancia de las variables usando la funcion importance()
## IncNodePurity
## CompPrice 196.79119
## Income 198.82432
## Advertising 236.38564
## Population 147.20301
## Price 540.66745
## ShelveLoc 538.52677
## Age 234.85851
## Education 105.39935
## Urban 21.58035
## US 32.67248
Para tener una visualización de dicha importancia, consideremos la siguiente grafica.
Note entonces que, al igual que sucedio con el bagging, las dos variables mas importantes para este modelo son ShelveLoc y Price, y su importancia se resalta de manera significativa sobre las demás variables.
A continuación, veamos el efecto de \(m\) sobre el bosque aleatorio. Crearemos una función que reciba un valor de \(m\) y calcule el error de validación.
get_mse_random_forest <- function(m, train, test) {
random_forest_test <- randomForest(Sales ~ ., data = train, mtry = m)
result_test <- predict(random_forest_test, newdata = test)
error <- mean((result_test - test$Sales) ^ 2)
return(error)
} Ahora graficaremos el valor del \(MSE\) con respecto a \(m\).
m_values <- 1:(ncol(carseats) - 1)
mse_values <- sapply(m_values, get_mse_random_forest, carseats_train, carseats_test)
m_mse_values = data.frame(m_values, mse_values)
g <- ggplot(data = m_mse_values) +
geom_line(mapping = aes(x = m_values, y = mse_values)) +
geom_point(mapping = aes(x = m_values, y = mse_values)) +
labs(x = "m", y = "MSE", title = "Valor de MSE para cada valor de m")
ggplotly(g)Luego, extraeremos el mejor \(m\) como sigue.
best_m <- m_mse_values[which(m_mse_values$mse_values == min(m_mse_values$mse_values)), "m_values"]
best_mse <- min(m_mse_values$mse_values)
best_m## [1] 7
De esta manera, se muestra que considerando 7 variables dentro de la decisión en cada nodo, se obtiene un \(MSE\) de 2.394646.
Sección 8.4 - Ejercicio 9
El conjunto de datos OJ contiene 1070 compras donde los clientes compraron un jugo de naranja de Citrus Hill o de Maid Orange Juice y se guardaron un número de características de los clientes.
a) Conjuntos de entrenamiento y validación
Crear un conjunto de entrenamiento con 800 observaciones y 270 de prueba
oj_train <- as_tibble(OJ) %>%
rownames_to_column() %>%
sample_n(size = 800)
oj_test <- as_tibble(OJ) %>%
rownames_to_column() %>%
anti_join(oj_train, by = "rowname")
oj_train <- oj_train %>% dplyr::select(-rowname)
oj_test <- oj_test %>% dplyr::select(-rowname)
dim(oj_train)## [1] 800 18
## [1] 270 18
b) Árbol de decisión
##
## Classification tree:
## tree(formula = Purchase ~ ., data = oj_train)
## Variables actually used in tree construction:
## [1] "LoyalCH" "PriceDiff" "WeekofPurchase" "ListPriceDiff"
## Number of terminal nodes: 8
## Residual mean deviance: 0.821 = 650.3 / 792
## Misclassification error rate: 0.1925 = 154 / 800
El árbol resultante tiene 10 nodos terminales, la variable más importante es LoyalCH con una ponderación de 63, con una tasa de clasificación incorrecta de 14.88% en entrenamiento.
c) Interepretar el resultado de uno de los nodos
Interpretemos el nodo 1: En este nodo ingresan 800 observaciones del conjunto de entrenamiento, hay 311 observaciones de la clase CH y 489 de la clase MM, se realiza una partición a la izquierda con 512 observaciones y la derecha con 285 observaciones bajo el criterio \(LoyalCH < 0.48285\).
d) Gráfica del árbol y resultados
Los nodos utilizan las variables CH y PriceDiff para realizar la clasificación, los nodos terminales indican la clasificación final.
e) Predicciones, matriz de confusión y clasificación correcta
oj_tree_class <- tibble(pred = predict(oj_tree, type="class", newdata = oj_test),
real = oj_test$Purchase)
oj_tree_class %>%
conf_mat(truth = real, estimate = pred) %>%
autoplot(type = "heatmap") +
labs(title = "Matriz de confusión para el modelo CART para clasificación")## # A tibble: 2 x 3
## .metric .estimator .estimate
## <chr> <chr> <dbl>
## 1 accuracy binary 0.830
## 2 kap binary 0.649
Se obtiene una tasa de clasificación correcta del 80%. En general la mayor cantidad de errores se cometió al clasificar las variables CH como MM, esto ocurrió 29 veces.
f) Aplicar cv.tree() para determinar el tamaño óptimo del árbol
## $size
## [1] 8 7 4 2 1
##
## $dev
## [1] 188 187 187 189 308
##
## $k
## [1] -Inf 0.0000000 0.6666667 2.5000000 146.0000000
##
## $method
## [1] "misclass"
##
## attr(,"class")
## [1] "prune" "tree.sequence"
g) Gráfica con el tamaño del árbol y la validación cruzada
h) Árbol con la menor tasa de error de validación cruzada incorreta
De acuerdo a la gráfica anterior, el árbol con tamaño 5 parece obtener el menor error de validación cruzada. Realicemos la poda de acuerdo a esta inforamción y grafiquemos el nuevo árbol.
i) Creación del árbol podado
j) Comparación de los errores de clasificación en entrenamiento
oj_prune_class_train <- tibble(pred = predict(oj_prune, type="class"),
real = oj_train$Purchase)
oj_prune_class_train %>%
conf_mat(truth = real, estimate = pred) %>%
autoplot(type = "heatmap") +
labs(title = "Matriz de confusión en entrenamiento para el modelo CART podado para clasificación")## # A tibble: 2 x 3
## .metric .estimator .estimate
## <chr> <chr> <dbl>
## 1 accuracy binary 0.808
## 2 kap binary 0.598
Se obtiene una clasificación correcta del 84.5% para el modelo CART podado, esto se compara con el modelo completo que da una clasificación correcta del 85.2%, es decir, el modelo con validación cruzada es ligeramente peor en el conjunto de entrenamiento pero tiene una mayor capacidad de generalización gracias a que evita sobreajustar los datos.
oj_tree_class_train <- tibble(pred = predict(oj_tree, type="class"),
real = oj_train$Purchase)
oj_tree_class_train %>%
metrics(truth = real, estimate = pred)## # A tibble: 2 x 3
## .metric .estimator .estimate
## <chr> <chr> <dbl>
## 1 accuracy binary 0.808
## 2 kap binary 0.598
k) Comparación los resultados del árbol podado vs el completo en validación
oj_prune_class <- tibble(pred = predict(oj_prune, type="class", newdata = oj_test),
real = oj_test$Purchase)
oj_prune_class %>%
conf_mat(truth = real, estimate = pred) %>%
autoplot(type = "heatmap") +
labs(title = "Matriz de confusión para el modelo CART podado para clasificación")## # A tibble: 2 x 3
## .metric .estimator .estimate
## <chr> <chr> <dbl>
## 1 accuracy binary 0.830
## 2 kap binary 0.649
En el árbol sin podar obtuvimos una tasa de clasificación correcta del 80%, sin embargo, al realizar la poda se obtuvo una tasa de clasificación correcta del 80.4%, que podría ser una mejora apreciable dependiendo del volumen de ventas de jugo de naranja asociado que trae la predicción correcta.
Sección 8.4 - Ejercicio 10
El objetivo de esta sección es usar boosting para hacer una predicción sobre la variable Salary en el conjunto de datos Hitters
El conjunto de datos recoge información sobre jugadores de béisbol de las ligas mayores de Estados Unidos. A continuación se muestra un resumen de dicho conjunto de datos.
## # A tibble: 6 x 20
## AtBat Hits HmRun Runs RBI Walks Years CAtBat CHits CHmRun CRuns CRBI
## <int> <int> <int> <int> <int> <int> <int> <int> <int> <int> <int> <int>
## 1 293 66 1 30 29 14 1 293 66 1 30 29
## 2 315 81 7 24 38 39 14 3449 835 69 321 414
## 3 479 130 18 66 72 76 3 1624 457 63 224 266
## 4 496 141 20 65 78 37 11 5628 1575 225 828 838
## 5 321 87 10 39 42 30 2 396 101 12 48 46
## 6 594 169 4 74 51 35 11 4408 1133 19 501 336
## # ... with 8 more variables: CWalks <int>, League <fct>, Division <fct>,
## # PutOuts <int>, Assists <int>, Errors <int>, Salary <dbl>, NewLeague <fct>
a) Remoción de NA’s y escalamiento logaritmico
Inicialmente se eliminarán aquellas filas donde el salario tenga el valor NA. Para ello primero observemos si la columna Salary es la única que tiene NA.
## AtBat Hits HmRun Runs
## Min. : 16.0 Min. : 1 Min. : 0.00 Min. : 0.00
## 1st Qu.:255.2 1st Qu.: 64 1st Qu.: 4.00 1st Qu.: 30.25
## Median :379.5 Median : 96 Median : 8.00 Median : 48.00
## Mean :380.9 Mean :101 Mean :10.77 Mean : 50.91
## 3rd Qu.:512.0 3rd Qu.:137 3rd Qu.:16.00 3rd Qu.: 69.00
## Max. :687.0 Max. :238 Max. :40.00 Max. :130.00
##
## RBI Walks Years CAtBat
## Min. : 0.00 Min. : 0.00 Min. : 1.000 Min. : 19.0
## 1st Qu.: 28.00 1st Qu.: 22.00 1st Qu.: 4.000 1st Qu.: 816.8
## Median : 44.00 Median : 35.00 Median : 6.000 Median : 1928.0
## Mean : 48.03 Mean : 38.74 Mean : 7.444 Mean : 2648.7
## 3rd Qu.: 64.75 3rd Qu.: 53.00 3rd Qu.:11.000 3rd Qu.: 3924.2
## Max. :121.00 Max. :105.00 Max. :24.000 Max. :14053.0
##
## CHits CHmRun CRuns CRBI
## Min. : 4.0 Min. : 0.00 Min. : 1.0 Min. : 0.00
## 1st Qu.: 209.0 1st Qu.: 14.00 1st Qu.: 100.2 1st Qu.: 88.75
## Median : 508.0 Median : 37.50 Median : 247.0 Median : 220.50
## Mean : 717.6 Mean : 69.49 Mean : 358.8 Mean : 330.12
## 3rd Qu.:1059.2 3rd Qu.: 90.00 3rd Qu.: 526.2 3rd Qu.: 426.25
## Max. :4256.0 Max. :548.00 Max. :2165.0 Max. :1659.00
##
## CWalks League Division PutOuts Assists
## Min. : 0.00 A:175 E:157 Min. : 0.0 Min. : 0.0
## 1st Qu.: 67.25 N:147 W:165 1st Qu.: 109.2 1st Qu.: 7.0
## Median : 170.50 Median : 212.0 Median : 39.5
## Mean : 260.24 Mean : 288.9 Mean :106.9
## 3rd Qu.: 339.25 3rd Qu.: 325.0 3rd Qu.:166.0
## Max. :1566.00 Max. :1378.0 Max. :492.0
##
## Errors Salary NewLeague
## Min. : 0.00 Min. : 67.5 A:176
## 1st Qu.: 3.00 1st Qu.: 190.0 N:146
## Median : 6.00 Median : 425.0
## Mean : 8.04 Mean : 535.9
## 3rd Qu.:11.00 3rd Qu.: 750.0
## Max. :32.00 Max. :2460.0
## NA's :59
Note entonces que la única columna que tiene NA es Salary. Ahora podemos eliminar estas filas y verificar que no quedan NA
## AtBat Hits HmRun Runs
## Min. : 19.0 Min. : 1.0 Min. : 0.00 Min. : 0.00
## 1st Qu.:282.5 1st Qu.: 71.5 1st Qu.: 5.00 1st Qu.: 33.50
## Median :413.0 Median :103.0 Median : 9.00 Median : 52.00
## Mean :403.6 Mean :107.8 Mean :11.62 Mean : 54.75
## 3rd Qu.:526.0 3rd Qu.:141.5 3rd Qu.:18.00 3rd Qu.: 73.00
## Max. :687.0 Max. :238.0 Max. :40.00 Max. :130.00
## RBI Walks Years CAtBat
## Min. : 0.00 Min. : 0.00 Min. : 1.000 Min. : 19.0
## 1st Qu.: 30.00 1st Qu.: 23.00 1st Qu.: 4.000 1st Qu.: 842.5
## Median : 47.00 Median : 37.00 Median : 6.000 Median : 1931.0
## Mean : 51.49 Mean : 41.11 Mean : 7.312 Mean : 2657.5
## 3rd Qu.: 71.00 3rd Qu.: 57.00 3rd Qu.:10.000 3rd Qu.: 3890.5
## Max. :121.00 Max. :105.00 Max. :24.000 Max. :14053.0
## CHits CHmRun CRuns CRBI
## Min. : 4.0 Min. : 0.00 Min. : 2.0 Min. : 3.0
## 1st Qu.: 212.0 1st Qu.: 15.00 1st Qu.: 105.5 1st Qu.: 95.0
## Median : 516.0 Median : 40.00 Median : 250.0 Median : 230.0
## Mean : 722.2 Mean : 69.24 Mean : 361.2 Mean : 330.4
## 3rd Qu.:1054.0 3rd Qu.: 92.50 3rd Qu.: 497.5 3rd Qu.: 424.5
## Max. :4256.0 Max. :548.00 Max. :2165.0 Max. :1659.0
## CWalks League Division PutOuts Assists
## Min. : 1.0 A:139 E:129 Min. : 0.0 Min. : 0.0
## 1st Qu.: 71.0 N:124 W:134 1st Qu.: 113.5 1st Qu.: 8.0
## Median : 174.0 Median : 224.0 Median : 45.0
## Mean : 260.3 Mean : 290.7 Mean :118.8
## 3rd Qu.: 328.5 3rd Qu.: 322.5 3rd Qu.:192.0
## Max. :1566.0 Max. :1377.0 Max. :492.0
## Errors Salary NewLeague
## Min. : 0.000 Min. : 67.5 A:141
## 1st Qu.: 3.000 1st Qu.: 190.0 N:122
## Median : 7.000 Median : 425.0
## Mean : 8.593 Mean : 535.9
## 3rd Qu.:13.000 3rd Qu.: 750.0
## Max. :32.000 Max. :2460.0
Ahora realicemos un escalamiento logaritmico sobre la columna Salary
## # A tibble: 6 x 20
## AtBat Hits HmRun Runs RBI Walks Years CAtBat CHits CHmRun CRuns CRBI
## <int> <int> <int> <int> <int> <int> <int> <int> <int> <int> <int> <int>
## 1 315 81 7 24 38 39 14 3449 835 69 321 414
## 2 479 130 18 66 72 76 3 1624 457 63 224 266
## 3 496 141 20 65 78 37 11 5628 1575 225 828 838
## 4 321 87 10 39 42 30 2 396 101 12 48 46
## 5 594 169 4 74 51 35 11 4408 1133 19 501 336
## 6 185 37 1 23 8 21 2 214 42 1 30 9
## # ... with 8 more variables: CWalks <int>, League <fct>, Division <fct>,
## # PutOuts <int>, Assists <int>, Errors <int>, Salary <dbl>, NewLeague <fct>
b) Seleccion de los conjuntos de prueba y validación
Se seleccionan las primeras 200 filas como el conjunto de prueba y las filas restantes se seleccionan como conjunto de validación, obteniendo los dos siguientes conjuntos de datos de manera respectiva.
## # A tibble: 200 x 20
## AtBat Hits HmRun Runs RBI Walks Years CAtBat CHits CHmRun CRuns CRBI
## <int> <int> <int> <int> <int> <int> <int> <int> <int> <int> <int> <int>
## 1 315 81 7 24 38 39 14 3449 835 69 321 414
## 2 479 130 18 66 72 76 3 1624 457 63 224 266
## 3 496 141 20 65 78 37 11 5628 1575 225 828 838
## 4 321 87 10 39 42 30 2 396 101 12 48 46
## 5 594 169 4 74 51 35 11 4408 1133 19 501 336
## 6 185 37 1 23 8 21 2 214 42 1 30 9
## 7 298 73 0 24 24 7 3 509 108 0 41 37
## 8 323 81 6 26 32 8 2 341 86 6 32 34
## 9 401 92 17 49 66 65 13 5206 1332 253 784 890
## 10 574 159 21 107 75 59 10 4631 1300 90 702 504
## # ... with 190 more rows, and 8 more variables: CWalks <int>, League <fct>,
## # Division <fct>, PutOuts <int>, Assists <int>, Errors <int>, Salary <dbl>,
## # NewLeague <fct>
## # A tibble: 63 x 20
## AtBat Hits HmRun Runs RBI Walks Years CAtBat CHits CHmRun CRuns CRBI
## <int> <int> <int> <int> <int> <int> <int> <int> <int> <int> <int> <int>
## 1 419 101 18 65 58 92 20 9528 2510 548 1509 1659
## 2 376 82 21 42 60 35 5 1770 408 115 238 299
## 3 486 145 11 51 76 40 11 3967 1102 67 410 497
## 4 246 76 5 35 39 13 6 912 234 12 102 96
## 5 205 52 8 31 27 17 12 5134 1323 56 643 445
## 6 348 90 11 50 45 43 10 2288 614 43 295 273
## 7 523 135 8 52 44 52 9 3368 895 39 377 284
## 8 312 68 2 32 22 24 1 312 68 2 32 22
## 9 496 119 8 57 33 21 7 3358 882 36 365 280
## 10 126 27 3 8 10 5 4 239 49 3 16 13
## # ... with 53 more rows, and 8 more variables: CWalks <int>, League <fct>,
## # Division <fct>, PutOuts <int>, Assists <int>, Errors <int>, Salary <dbl>,
## # NewLeague <fct>
c) Aplicación de boosting y determinación del parametro \(\lambda\) de contracción sobre el conjunto de entrenamiento.
En esta sección realizaremos boosting y calcularemos el mejor parámetro de contracción \(\lambda\) sobre un determinado rango. Para esto, crearemos una función que calcule el \(MSE\) sobre el conjunto de entrenamiento recibiendo \(\lambda\) como parámetro.
get_mse_boosting_train <- function(lambda, train) {
test_gbm <- gbm(Salary ~ ., data = train, distribution = "gaussian", n.trees = 1000, shrinkage = lambda)
predict_train <- predict(test_gbm, newdata = train, n.trees = 1000)
error <- mean((predict_train - train$Salary) ^ 2)
return(error)
}Ahora definamos el rango sobre el cual se quiere graficar los \(MSE\), y realicemos dicha gráfica usando la función que se definió anteriormente.
lambda_values = seq(0.01, 1, by = 0.01)
mse_train_values_hitters <- sapply(lambda_values, get_mse_boosting_train, hitters_train)
lambda_mse_train_values_df <- data.frame(lambda_values, mse_values = mse_train_values_hitters)
best_lambda_train_boosting <- lambda_mse_train_values_df[which(lambda_mse_train_values_df$mse_values == min(lambda_mse_train_values_df$mse_values)), "lambda_values"]
best_mse_train_boosting <- min(lambda_mse_train_values_df$mse_values)
ggplot(data = lambda_mse_train_values_df) +
geom_line(mapping = aes(x = lambda_values, y = mse_values)) +
labs(x = TeX("$\\lambda$"), y = "MSE de entrenamiento", title = "MSE de entrenamiento para diferentes valores de contracción") +
theme_economist()En esta gráfica se puede apreciar que el mejor valor para \(\lambda\) es 0.98 con el cual se obtiene un \(MSE\) de 1.15062610^{-4}.
d) Determinación del parametro \(\lambda\) de contracción sobre el conjunto de prueba
Ahora graficaremos el \(MSE\) de validación para diferentes valores de \(\lambda\) y analizaremos cual de ellos obtiene un valor mas apropiado.
Para ello, primero debemos construir una función que permita calcular el \(MSE\) de validación para un determinado valor de \(\lambda\),
get_mse_boosting_test <- function(lambda, train, test) {
test_gbm <- gbm(Salary ~ ., data = train, distribution = "gaussian", n.trees = 1000, shrinkage = lambda)
predict_test <- predict(test_gbm, newdata = test, n.trees = 1000)
error <- mean((predict_test - test$Salary) ^ 2)
return(error)
}Ahora grafiquemos el correspondiente valor del \(MSE\) para cada uno de los posibles valores de \(\lambda\) seleccionados.
mse_test_values_hitters <- sapply(lambda_values, get_mse_boosting_test, hitters_train, hitters_test)
lambda_test_mse_values_df <- data.frame(lambda_values, mse_values = mse_test_values_hitters)
best_lambda_test_boosting <- lambda_test_mse_values_df[which(lambda_test_mse_values_df$mse_values == min(lambda_test_mse_values_df$mse_values)), "lambda_values"]
best_mse_test_boosting <- min(lambda_test_mse_values_df$mse_values)
ggplot(data = lambda_test_mse_values_df) +
geom_line(mapping = aes(x = lambda_values, y = mse_values)) +
labs(x = TeX("$\\lambda$"), y = "MSE de validación", title = "MSE de validación para diferentes valores de contracción") +
theme_economist()Note entonces que el mejor valor para \(\lambda\) es 0.08, del cual se obtiene un valor para el \(MSE\) de validación de 0.0480158.
Como conclusión respecto a los valore encontrados en la gráfica del literal anterior, se puede observar entonces que al aumentar el \(\lambda\), el \(MSE\) de entrenamiento desciende ya que el algoritmo de optimización pretende minimizar el \(MSE\) de los datos con el cuál se entrena el modelo. Pero algo muy interesante que hay que notar, es que contrario a lo que pasa en el conjunto de entrenamiento, el \(MSE\) de validación aumenta a medida que aumentamos \(\lambda\), lo que sugiere que hay un sobreajuste del modelo para valores de \(\lambda\) cercanos a 1 o mayores. Esto muestra que juzgar un modelo mediante un conjunto de entrenamiento no es adecuado, pues pueden haber parámetros que no esten generalizando bien el comportamiento real de los datos.
e) Comparación de boosting con métodos previos
En esta sección compararemos el \(MSE\) de boosting respecto a otros métodos vistos ateriormente como regresión lineal múltiple y regresión de componentes principales.
Primero calculemos el \(MSE\) para el metodo de regresión lineal múltiple
linear_hitters <- lm(Salary ~ ., data = hitters_train)
linear_hitters_pred <- predict(linear_hitters, hitters_test)
error_linear <- mean((linear_hitters_pred - hitters_test$Salary))
error_linear## [1] 0.0188556
Por otro lado, apliquemos regresión con componentes principales para extraer el \(MSE\).
pcr_hitters <- pcr(Salary ~ ., data = hitters_train, scale = TRUE, validation = "CV")
summary(pcr_hitters)## Data: X dimension: 200 19
## Y dimension: 200 1
## Fit method: svdpc
## Number of components considered: 19
##
## VALIDATION: RMSEP
## Cross-validated using 10 random segments.
## (Intercept) 1 comps 2 comps 3 comps 4 comps 5 comps 6 comps
## CV 0.3981 0.2796 0.2808 0.2810 0.2809 0.2797 0.2787
## adjCV 0.3981 0.2792 0.2803 0.2805 0.2804 0.2788 0.2780
## 7 comps 8 comps 9 comps 10 comps 11 comps 12 comps 13 comps
## CV 0.2769 0.2768 0.2813 0.2830 0.2854 0.2834 0.2807
## adjCV 0.2761 0.2759 0.2802 0.2818 0.2841 0.2821 0.2796
## 14 comps 15 comps 16 comps 17 comps 18 comps 19 comps
## CV 0.2811 0.2816 0.2787 0.2777 0.2795 0.2800
## adjCV 0.2797 0.2801 0.2771 0.2760 0.2777 0.2782
##
## TRAINING: % variance explained
## 1 comps 2 comps 3 comps 4 comps 5 comps 6 comps 7 comps 8 comps
## X 39.40 60.67 70.76 79.64 84.73 89.29 92.38 95.02
## Salary 52.42 52.67 53.22 53.56 54.56 55.34 55.70 56.48
## 9 comps 10 comps 11 comps 12 comps 13 comps 14 comps 15 comps
## X 96.32 97.26 98.06 98.72 99.23 99.54 99.77
## Salary 56.68 56.77 56.80 57.41 57.91 58.68 58.75
## 16 comps 17 comps 18 comps 19 comps
## X 99.90 99.97 99.99 100.00
## Salary 60.27 61.27 61.27 61.49
En este caso, encontramos que el valor óptimo se logró con \(M = 17\) por tanto calcularemos el error de validación con dicho número de componentes, obteniendo así el siguiente \(MSE\).
pcr_hitters_pred <- predict(pcr_hitters, hitters_test, ncomp = 17)
error_pcr <- mean((pcr_hitters_pred - hitters_test$Salary))
error_pcr## [1] 0.02113112
Finalmente, concluimos que en este caso el boosting no era la mejor opción, pues encontramos que el \(MSE\) para la regresión lineal multiple era mucho menor, y de hecho es el menor de los tres \(MSE\) de prueba calculados.
f) Importancia de las variables
Ahora observaremos cual es la improtancia de las variables respecto a este método. Para ello entrenemos el modelo con el mejor \(\lambda\) obtenido y analicemos la tabla de importancias.
gbm_hitters <- gbm(Salary ~ ., data = hitters_train, distribution = "gaussian", n.trees = 1000, shrinkage = best_mse_test_boosting)
summary(gbm_hitters)## var rel.inf
## CAtBat CAtBat 24.3085290
## CWalks CWalks 8.8695835
## CRBI CRBI 6.8122643
## Walks Walks 6.6965360
## Years Years 6.5480806
## CHits CHits 5.8999306
## PutOuts PutOuts 5.8424749
## CHmRun CHmRun 5.3895812
## RBI RBI 5.3256242
## CRuns CRuns 5.2700070
## Hits Hits 3.5469571
## Assists Assists 3.4812787
## HmRun HmRun 3.0232305
## Errors Errors 2.6757639
## Runs Runs 2.4816016
## AtBat AtBat 2.4331138
## Division Division 0.7205061
## League League 0.3461861
## NewLeague NewLeague 0.3287510
Se puede concluir entonces que las variables CAtBat, CHits y CRuns son las variables que más influyentes dentro del modelo.
g) Aplicación de bagging
En este caso aplicaremos bagging sobre este conjunto de datos y observaremos su \(MSE\) de prueba. Dicho \(MSE\) se muestra a continuación
bagging_hitters <- randomForest(Salary ~ ., data = hitters_train, mtry = ncol(hitters_train) - 1, importance = TRUE)
bagging_predict_test_hitters <- predict(bagging_hitters, hitters_test)
bagging_error_mse <- mean((bagging_predict_test_hitters - hitters_test$Salary) ^ 2)
bagging_error_mse## [1] 0.04338928
Como conclusión podemos observar que el bagging no ofrece un cambio positivo respecto al error, pues al analizar los \(MSE\) de prueba obtenidos con los anteriores métodos, el método lineal sigue siendo el más efectivo de todos.
Sección 8.4 - Ejercicio 11
El conjunto de datos Caravan contiene 5822 registros de clientes de compañía de seguros, donde cada uno de estos se posee 86 variables con datos sociodemográficos y sus productos asociados.
La variable Purchase indica si el cliente compró o no un programa de seguros de la compañía.
a) Conjunto de entrenamiento con las primeras 1000 observaciones
Dado que el problema a resolver es de clasificación, se utiliza el argumento distribution = "bernoulli", sin embargo esta respuesta requiere que transformemos la variable respuesta de “Yes” y “No” a 1 y 0.
b) Ajuste de un modelo Random Forest con Boosting utilizando \(\lambda = 0.01\), 1000 árboles y variables más importantes
caravan_boost <- gbm(Purchase ~ ., data = caravan_train, distribution = "bernoulli", n.trees = 1000, shrinkage = 0.01)
head(summary(caravan_boost), n = 10)## var rel.inf
## PPERSAUT PPERSAUT 15.137917
## MKOOPKLA MKOOPKLA 9.882042
## MOPLHOOG MOPLHOOG 7.600803
## MINK3045 MINK3045 5.043398
## MBERMIDD MBERMIDD 4.814501
## PBRAND PBRAND 4.770935
## MGODGE MGODGE 4.521611
## ABRAND ABRAND 4.268103
## MOSTYPE MOSTYPE 3.144056
## PWAPART PWAPART 2.808144
Las variables más importantes del modelo son PPERSAUT y MKOOPKLA, lastimosamente no hay una descripción clara de qué representan en el conjunto de datos.
En las gráficas de importancia parcial para estas dos variables se ilustra su efecto marginal sobre la respuesta luego de integrar las demás variables. En este caso se identifica que la probabilidad de que ocurra una compra decrece con los valores de ambas variables.
c) Predicción en el conjunto de prueba para probabilidades mayores a 20%
caravan_boost_class <- tibble(value = predict(caravan_boost, newdata = caravan_test, n.trees = 1000, type = "response"),
pred = as_factor(case_when(value > 0.2 ~ 1,
TRUE ~ 0)),
real = as_factor(caravan_test$Purchase))
caravan_boost_class %>%
conf_mat(truth = real, estimate = pred) %>%
autoplot(type = "heatmap") +
labs(title = "Matriz de confusión para el modelo de Random Forest con Boosting")## # A tibble: 2 x 3
## .metric .estimator .estimate
## <chr> <chr> <dbl>
## 1 accuracy binary 0.921
## 2 kap binary 0.113
Se obtiene una tasa de clasificación correcta del 92%, sin embargo, sólo el 11.7% de los clientes que sí realizaron la compra se predicen correctamente, lo cual puede ser un resultado inadecuado dependiendo de las decisiones de negocio realizadas con esta predicción. Esto se debe a que la clase “No” tiene muchos más valores y por ende la tasa de clasificación correcta no es una medida representativa para la clase “Yes”.
Al realizar comparación con la regresión logística se obtienen los siguientes resultados:
caravan_logistic <- glm(Purchase ~ ., data = caravan_train, family = binomial)
caravan_logistic_class <- tibble(value = predict(caravan_logistic, newdata = caravan_test, type = "response"),
pred = as_factor(case_when(value > 0.2 ~ 1,
TRUE ~ 0)),
real = as_factor(caravan_test$Purchase))
caravan_logistic_class %>%
conf_mat(truth = real, estimate = pred) %>%
autoplot(type = "heatmap") +
labs(title = "Matriz de confusión para el modelo de regresión logística")## # A tibble: 2 x 3
## .metric .estimator .estimate
## <chr> <chr> <dbl>
## 1 accuracy binary 0.880
## 2 kap binary 0.104
La tasa de clasificación correcta del modelo de regresión logística disminuye al 88% de 92%, sin embargo, ésta discrimina mejor los casos en los que sí se realiza una compra, pasando del 11.7% del modelo de boosting a un 20%.
Sección 8.4 - Ejercicio 12
En este ejercicio utilizaremos información médica de pacientes para predecir si tienen o no enfermedad del corazón. El objetivo es predecir la variable hd, donde 0 es saludable y 1 significa que padece una enfermedad del corazón.
heart_disease <- "http://archive.ics.uci.edu/ml/machine-learning-databases/heart-disease/processed.cleveland.data" %>%
read_csv(col_names = FALSE) ## Parsed with column specification:
## cols(
## X1 = col_double(),
## X2 = col_double(),
## X3 = col_double(),
## X4 = col_double(),
## X5 = col_double(),
## X6 = col_double(),
## X7 = col_double(),
## X8 = col_double(),
## X9 = col_double(),
## X10 = col_double(),
## X11 = col_double(),
## X12 = col_character(),
## X13 = col_character(),
## X14 = col_double()
## )
names(heart_disease) = c(
"age",
"sex",# 0 = female, 1 = male
"cp", # chest pain
# 1 = typical angina,
# 2 = atypical angina,
# 3 = non-anginal pain,
# 4 = asymptomatic
"trestbps", # resting blood pressure (in mm Hg)
"chol", # serum cholestoral in mg/dl
"fbs", # fasting blood sugar if less than 120 mg/dl, 1 = TRUE, 0 = FALSE
"restecg", # resting electrocardiographic results
# 1 = normal
# 2 = having ST-T wave abnormality
# 3 = showing probable or definite left ventricular hypertrophy
"thalach", # maximum heart rate achieved
"exang", # exercise induced angina, 1 = yes, 0 = no
"oldpeak", # ST depression induced by exercise relative to rest
"slope", # the slope of the peak exercise ST segment
# 1 = upsloping
# 2 = flat
# 3 = downsloping
"ca", # number of major vessels (0-3) colored by fluoroscopy
"thal", # this is short of thalium heart scan
# 3 = normal (no cold spots)
# 6 = fixed defect (cold spots during rest and exercise)
# 7 = reversible defect (when cold spots only appear during exercise)
"hd" # (the predicted attribute) - diagnosis of heart disease
# 0 if less than or equal to 50% diameter narrowing
# 1 if greater than 50% diameter narrowing
)
heart_disease <- heart_disease %>%
mutate(
sex = case_when(sex == 1 ~ "Male",
sex == 0 ~ "Female"),
sex = as_factor(sex),
cp = as_factor(cp),
fbs = as_factor(fbs),
restecg = as_factor(restecg),
exang = as_factor(exang),
slope = as_factor(slope),
ca = as.integer(ca) %>% as_factor(),
thal = as.integer(thal) %>% as_factor(),
hd = case_when(hd == 0 ~ 0,
hd >= 1 ~ 1)) %>%
drop_na()## Warning in eval(lhs, parent, parent): NAs introducidos por coerción
## Warning in eval(lhs, parent, parent): NAs introducidos por coerción
## # A tibble: 6 x 14
## age sex cp trestbps chol fbs restecg thalach exang oldpeak slope
## <dbl> <fct> <fct> <dbl> <dbl> <fct> <fct> <dbl> <fct> <dbl> <fct>
## 1 63 Male 1 145 233 1 2 150 0 2.3 3
## 2 67 Male 4 160 286 0 2 108 1 1.5 2
## 3 67 Male 4 120 229 0 2 129 1 2.6 2
## 4 37 Male 3 130 250 0 0 187 0 3.5 3
## 5 41 Fema~ 2 130 204 0 2 172 0 1.4 1
## 6 56 Male 2 120 236 0 0 178 0 0.8 1
## # ... with 3 more variables: ca <fct>, thal <fct>, hd <dbl>
hd_split_ratio <- 0.75
hd_smp_size <- floor(hd_split_ratio * nrow(heart_disease))
hd_train_ind <- sample(seq_len(nrow(heart_disease)), size = hd_smp_size)
hd_train <- heart_disease[hd_train_ind, ]
hd_test <- heart_disease[-hd_train_ind, ]Modelo de Boosting para clasificación de enfermedad del corazón
Las 3 variables que más influyen en el diagnóstico correcto de la enfermedad del corazón son ca, que es el número de major vessels por la técnica de fluoroscopía, thal, que es el escaneo de talio del corazón y cp, que corresponde al tipo de dolor de pecho.
heart_boost <- gbm(hd ~ ., data = hd_train, distribution = "bernoulli", n.trees = 1000, shrinkage = 0.01)
heart_boost## gbm(formula = hd ~ ., distribution = "bernoulli", data = hd_train,
## n.trees = 1000, shrinkage = 0.01)
## A gradient boosted model with bernoulli loss function.
## 1000 iterations were performed.
## There were 13 predictors of which 13 had non-zero influence.
## var rel.inf
## cp cp 21.92997175
## ca ca 20.70247146
## thal thal 12.63405880
## thalach thalach 12.52643882
## oldpeak oldpeak 10.09063206
## trestbps trestbps 4.69332333
## slope slope 4.68476146
## age age 4.50876171
## chol chol 4.15596369
## sex sex 2.45492073
## exang exang 1.25751066
## restecg restecg 0.32702078
## fbs fbs 0.03416475
heart_boost_class <- tibble(value = predict(heart_boost, newdata = hd_test, n.trees = 1000, type = "response"),
pred = as_factor(case_when(value > 0.9 ~ 1,
TRUE ~ 0)),
real = as_factor(hd_test$hd))
heart_boost_class %>%
conf_mat(truth = real, estimate = pred) %>%
autoplot(type = "heatmap") +
labs(title = "Matriz de confusión para el modelo de Boosting")## # A tibble: 2 x 3
## .metric .estimator .estimate
## <chr> <chr> <dbl>
## 1 accuracy binary 0.68
## 2 kap binary 0.313
El modelo de boosting clasificó de manera correcta el 74.7% de las observaciones del conjunto de prueba, con sólo 1 falso negativo. Es importante notar que en un caso médico como éste vale la pena crear falsos positivos, ya que son menos costosos, sin embargo, no diagnosticar a un paciente con enfermedad del corazón podría tener consecuencias fatales.
Por este motivo, se ha cambiado el umbral de probabilidad al 90% con el objetivo de producir menos falsos negativos, sacrificando parte de la capacidad de clasificación correcta de falsos positivos.
Modelo de Bagging para clasificación de enfermedad del corazón
Las 3 variables que más influyen en el diagnóstico correcto de la enfermedad del corazón son ca, que es el número de major vessels por la técnica de fluoroscopía, thal, que es el escaneo de talio del corazón y cp, que corresponde al tipo de dolor de pecho.
hd_train <- hd_train %>%
mutate(hd = as_factor(hd))
hd_test <- hd_test %>%
mutate(hd = as_factor(hd))
heart_bag <- randomForest(hd ~ ., data = hd_train, importance = TRUE, mtry = ncol(hd_train)-1)
heart_bag##
## Call:
## randomForest(formula = hd ~ ., data = hd_train, importance = TRUE, mtry = ncol(hd_train) - 1)
## Type of random forest: classification
## Number of trees: 500
## No. of variables tried at each split: 13
##
## OOB estimate of error rate: 18.92%
## Confusion matrix:
## 0 1 class.error
## 0 102 17 0.1428571
## 1 25 78 0.2427184
heart_bag_class <- tibble(pred = predict(heart_bag, newdata = hd_test),
real = as_factor(hd_test$hd))
heart_bag_class %>%
conf_mat(truth = real, estimate = pred) %>%
autoplot(type = "heatmap") +
labs(title = "Matriz de confusión para el modelo de Random Forest con Bagging")## # A tibble: 2 x 3
## .metric .estimator .estimate
## <chr> <chr> <dbl>
## 1 accuracy binary 0.813
## 2 kap binary 0.621
El modelo de bagging tiene un rendimiento considerablemente menor con una tasa de clasificación correcta del 77.3%, especialmente porque no puede ser parametrizado para minimizar los falsos negativos, lo cuál podría resultar en personas no diagnosticadas de manera peligrosa.
Modelo de Random Forest para clasificación de enfermedad del corazón
Las 3 variables que más influyen en el diagnóstico correcto de la enfermedad del corazón son ca, que es el número de major vessels por la técnica de fluoroscopía, thal, que es el escaneo de talio del corazón y cp, que corresponde al tipo de dolor de pecho.
##
## Call:
## randomForest(formula = hd ~ ., data = hd_train, importance = TRUE)
## Type of random forest: classification
## Number of trees: 500
## No. of variables tried at each split: 3
##
## OOB estimate of error rate: 18.47%
## Confusion matrix:
## 0 1 class.error
## 0 102 17 0.1428571
## 1 24 79 0.2330097
heart_rf_class <- tibble(pred = predict(heart_rf, newdata = hd_test),
real = as_factor(hd_test$hd))
heart_rf_class %>%
conf_mat(truth = real, estimate = pred) %>%
autoplot(type = "heatmap") +
labs(title = "Matriz de confusión para el modelo de Random Forest")## # A tibble: 2 x 3
## .metric .estimator .estimate
## <chr> <chr> <dbl>
## 1 accuracy binary 0.827
## 2 kap binary 0.648
El modelo de Random Forst sin Bagging tiene un rendimiento ligeramente mejor que el modelo con Bagging, sin embargo sufre del mismo problema de clasificación de falsos negativos altos, lo qeu podría conducir a fatalidades.
Modelo de regresión logística para clasificación de enfermedad del corazón
heart_log <- glm(hd ~ ., data = hd_train, family = "binomial")
heart_log_class <- tibble(value = predict(heart_log, newdata = hd_test, type = "response"),
pred = as_factor(case_when(value > 0.8 ~ 1,
TRUE ~ 0)),
real = as_factor(hd_test$hd))
heart_log_class %>%
conf_mat(truth = real, estimate = pred) %>%
autoplot(type = "heatmap") +
labs(title = "Matriz de confusión para el modelo de Regresión logística")## # A tibble: 2 x 3
## .metric .estimator .estimate
## <chr> <chr> <dbl>
## 1 accuracy binary 0.76
## 2 kap binary 0.498
En este caso, incrementar el umbral a 90% no mejora los falsos negativos, por lo que se deja en 80%, lo que produce una tasa de clasificación correcta del 85.3%, presentando mejor rendimiento que las técnicas de Boosting, Bagging y Random Forest.
Sección 9.7 - Ejercicio 4
El objetivo de este ejercicio es verificar la efectividad del metodo SVM para clasificar un conjunto de datos donde no hay una separación lineal entre ellos. Observaremos entonces una comparación entre un clasificador de soporte vectorial y una máquina de soporte vectorial con kernel polinomial o con un kernel radial.
Primero crearemos el conjunto de datos sintéticos con una separación no lineal.
x1_values <- runif(100, min = 0, max = 1)
x2_values <- runif(100, min = 0, max = 1)
radius <- 3 / 8
syntetic_df <- as_tibble(data.frame(x1 = x1_values, x2 = x2_values))
syntetic_df <- syntetic_df %>% mutate(
Class = case_when(
(x1 - 0.5) ^ 2 + (x2 - 0.5) ^ 2 <= radius ^ 2 ~ 1,
(x1 - 0.5) ^ 2 + (x2 - 0.5) ^ 2 > radius ^ 2 ~ 0,
),
Class = as.factor(Class)
)
head(syntetic_df)## # A tibble: 6 x 3
## x1 x2 Class
## <dbl> <dbl> <fct>
## 1 0.309 0.913 0
## 2 0.875 0.141 0
## 3 0.657 0.405 1
## 4 0.709 0.760 1
## 5 0.432 0.188 1
## 6 0.709 0.152 0
Grafiquemos el conjunto de datos generado.
ggplot(data = syntetic_df) +
geom_point(mapping = aes(x = x1, y = x2, color = Class)) +
geom_circle(mapping = aes(x0 = 0.5, y0 = 0.5, r = radius)) +
theme_economist() +
labs(x = TeX("$X_1$"), y = TeX("$X_2$"), title = "Datos sintéticos")Ahora dividiremos el conjunto de datos en un conjunto para entrenamiento y otro conjunto para validación, de la siguiente manera.
set.seed(10)
split_ratio <- 0.75
smp_size <- floor(split_ratio * nrow(syntetic_df))
train_ind <- sample(seq_len(nrow(syntetic_df)), size = smp_size)
syntetic_train <- syntetic_df[train_ind, ]
syntetic_test <- syntetic_df[-train_ind, ]
dim(syntetic_train)## [1] 75 3
## [1] 25 3
Primero entrenemos un SVC (Support Vector Classifier) y vamos a extraer su error en el conjunto de entrenamiento. Es importante notar que tomaremos el costo para la frontera de clasificación como 10 y observaremos el de este costo sobre los demás modelos.
linear_svm_syntethic <- svm(Class ~ ., data = syntetic_train, kernel = "linear", cost = 10, scale = FALSE)
summary(linear_svm_syntethic)##
## Call:
## svm(formula = Class ~ ., data = syntetic_train, kernel = "linear",
## cost = 10, scale = FALSE)
##
##
## Parameters:
## SVM-Type: C-classification
## SVM-Kernel: linear
## cost: 10
##
## Number of Support Vectors: 60
##
## ( 31 29 )
##
##
## Number of Classes: 2
##
## Levels:
## 0 1
Grafiquemos la frontera de decisión y observemos los resultados.
Y ahora, vamos a extraer el error sobre el conjunto de entrenamiento.
Note entonces que la region clasifica de manera correcta el 61.3333333% de los datos de entrenamiento (esta es la tasa de clasificación correcta de entrenamiento para dicho modelo).
A continuación, computemos la tasa de clasificación correcta para el conjunto de validación.
predict_svm_syntetic_test <- predict(linear_svm_syntethic, syntetic_test)
linear_svm_test_error <- mean(predict_svm_syntetic_test == syntetic_test$Class)Se puede notar entonces que para el conjunto de datos de validación, se tiene una tasa de clasificación correcta del 64%. Note que en este caso la grafica es aproximada a la de entrenamiento, y además este clasificador es muy similar a tomar cualquier pareja \(X_1\) y \(X_2\) y asignarles de manera aleatoria uniforme la etiqueta 1 o 0.
Esto se puede ver un poco mas claro observando la matriz de confusión para este modelo.
svm_syntetic_test_class <- tibble(pred = predict_svm_syntetic_test ,
real = syntetic_test$Class)
svm_syntetic_test_class %>% conf_mat(truth = real, estimate = pred) %>%
autoplot(type = "heatmap") +
labs(title = "Matriz de confusión para el modelo SVM con kernel lineal")Veamos el comportamiento de un kernel radial para este conjunto de datos. Para este caso, tomaremos un valor de \(\gamma = 1\).
radial_svm_syntethic <- svm(Class ~ ., data = syntetic_train, kernel = "radial", gamma = 1, cost = 10, scale = FALSE)
summary(radial_svm_syntethic)##
## Call:
## svm(formula = Class ~ ., data = syntetic_train, kernel = "radial",
## gamma = 1, cost = 10, scale = FALSE)
##
##
## Parameters:
## SVM-Type: C-classification
## SVM-Kernel: radial
## cost: 10
##
## Number of Support Vectors: 33
##
## ( 17 16 )
##
##
## Number of Classes: 2
##
## Levels:
## 0 1
Realicemos una gráfica para este clasificador sobre el conjunto de entrenamiento, como se muestra a continuación.
Es importante apreciar que el clasificador radial hace una clasificación casi perfecta, o al menos así se muestra gráficamente, además toma un conjunto mucho mas razonable de vectores de soporte, que son los que están más cercanos a la frontera de clasificación, al contrario del metodo lineal, el cual tomaba vectores de soporte que no corresponden a la intuición del modelo.
Ahora, para verificar el argumento que se mostró en el parrafo anterior, verifiquemos el error de entrenamiento de este modelo.
Note entonces que la tasa de clasificación correcta para el conjunto de entrenamiento tomando un kernel radial es de 93.3333333%, lo cual muestra una efectividad superior al clasificador de soporte vectorial (usando kernel lineal)
De igual manera, observemos la tasa de clasificación correcta para este modelo con kernel radial para el conjunto de validación.
predict_svm_syntetic_test_radial <- predict(radial_svm_syntethic, syntetic_test)
radial_svm_test_error <- mean(predict_svm_syntetic_test_radial == syntetic_test$Class)
radial_svm_test_error## [1] 0.84
Observemos como es la matriz de confusión para este modelo con kernel radial.
svm_syntetic_test_radial_class <- tibble(pred = predict_svm_syntetic_test_radial ,
real = syntetic_test$Class)
svm_syntetic_test_radial_class %>% conf_mat(truth = real, estimate = pred) %>%
autoplot(type = "heatmap") +
labs(title = "Matriz de confusión para el modelo SVM con kernel radial")Note entonces que en este caso, la tasa de clasificación correcta es de 84%, lo que muestra que el kernel radial es superior al kernel lineal cuando no hay una separación lineal entre las clases, como era de esperarse. Sin embargo, el cambio entre ambos modelos es muy significativo, pues se pasa de tener un modelo cuya selección en la práctica corresponde a un clasificador aleatorio uniforme, a tener un modelo con un desempeño casi perfecto.
Sección 9.7 - Ejercicio 5
En esta sección se observará que mediante la aplicación de regresión logística, también se puede obtener una frontera de decisión no lineal como se logra con las máquinas de soporte vectorial.
a) Construcción de los datos sintéticos
Inicialmente vamos a generar un conjunto de datos sintéticos, con una separación cuadrática entre las clases. A continuación mostraremos la generación de este conjunto de datos.
x1_values <- runif(n = 500) - 0.5
x2_values <- runif(n = 500) - 0.5
focus <- 1/8
y <- 1*(x1_values ^ 2 - x2_values ^ 2 > focus ^ 2)
syntetic_quadratic <- tibble(x1 = x1_values, x2 = x2_values, y = as.factor(y))
head(syntetic_quadratic)## # A tibble: 6 x 3
## x1 x2 y
## <dbl> <dbl> <fct>
## 1 -0.439 0.248 1
## 2 -0.131 0.223 0
## 3 -0.0924 0.0458 0
## 4 0.414 -0.236 1
## 5 0.389 -0.194 1
## 6 0.484 0.286 1
Ahora construiremos los conjuntos de entrenamiento y de pruebas para entrenar los posteriores modelos.
set.seed(8)
split_ratio <- 0.75
smp_size <- floor(split_ratio * nrow(syntetic_quadratic))
train_ind <- sample(seq_len(nrow(syntetic_quadratic)), size = smp_size)
syntetic_quadratic_train <- syntetic_quadratic[train_ind, ]
syntetic_quadratic_test <- syntetic_quadratic[-train_ind, ]
dim(syntetic_quadratic_train)## [1] 375 3
## [1] 125 3
b) Gráfica de los datos sintéticos
Ahora veamos una gráfica del conjunto de datos para visualizar la separación de los datos.
g <- ggplot(data = syntetic_quadratic) +
geom_point(mapping = aes(x = x1, y = x2, color = y)) +
labs(x = "X1", y = "X2", color = "y", title = "Conjunto de datos sintéticos") +
theme_economist()
ggplotly(g)Note entonces que los elementos que tienen el valor de \(y = 1\) son aquellos elementos que verifican la ecuación \(x_1^2 - x_2^2 > r^2\) donde \(r = 1/8\), formando así una frontera de decisión cuadrática.
c) Entrenamiento de un modelo logístico
Tomando en cuenta los datos construidos en los literales anteriores, vamos a entrenar un modelo de regresión logística, usando \(X_1\) y \(X_2\) como predictores para \(Y\).
set.seed(16)
glm_syntetic_quadratic <- glm(y ~ ., data = syntetic_quadratic_train, family = "binomial")
summary(glm_syntetic_quadratic)##
## Call:
## glm(formula = y ~ ., family = "binomial", data = syntetic_quadratic_train)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.3021 -1.1018 -0.9749 1.2265 1.4325
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.1822 0.1046 -1.741 0.0816 .
## x1 0.4784 0.3542 1.351 0.1768
## x2 -0.4844 0.3557 -1.362 0.1733
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 516.95 on 374 degrees of freedom
## Residual deviance: 513.32 on 372 degrees of freedom
## AIC: 519.32
##
## Number of Fisher Scoring iterations: 4
d) Gráficas de entrenamiento
Ahora apiquemos el modelo al conjunto de entrenamientoy grafiquemos la predicción de la clasificación para este conjunto de datos.
set.seed(16)
predict_train_syntetic_quadratic <- predict(glm_syntetic_quadratic, syntetic_quadratic_train, type = "response")
response_predict_train_syntetic_quadratic <- rep(0, nrow(syntetic_quadratic_train))
response_predict_train_syntetic_quadratic[predict_train_syntetic_quadratic > 0.5] = 1
ggplot(data = data.frame(x1 = syntetic_quadratic_train$x1,
x2 = syntetic_quadratic_train$x2,
y = as.factor(response_predict_train_syntetic_quadratic))) +
geom_point(mapping = aes(x = x1, y = x2, color = y)) +
theme_economist() +
labs(x = TeX("X_1"), y = TeX("X_2"), color = "y", title = "Predicciones sobre el conjunto de entrenamiento")Note entonces que la frontera de clasificación esta en una linea que toma solo algunos puntos de la clase 1 y el resto de puntos los toma como clase 0. Claramente, dentro del conjunto de entrenamiento, esto no es una buena clasificación para este conjunto de datos.
La calidad del modelo se verá reflejada en la matriz de confusión sobre el conjunto de evaluación, como se muestra en la siguiente figura.
predict_test_syntetic_quadratic <- predict(glm_syntetic_quadratic, syntetic_quadratic_test, type = "response")
response_predict_test_syntetic_quadratic <- rep(0, nrow(syntetic_quadratic_test))
response_predict_test_syntetic_quadratic[predict_test_syntetic_quadratic > 0.5] = 1
glm_test_syntetic_quadratic_class <- tibble(truth = as.factor(syntetic_quadratic_test$y),
estimate = as.factor(response_predict_test_syntetic_quadratic))
glm_test_syntetic_quadratic_class %>% conf_mat(truth = truth, estimate = estimate) %>%
autoplot(type = "heatmap") +
labs(title = "Matriz de confusión para el modelo logístico")En este punto, la matriz de confusión revela que el modelo tiene especial dificultad en el reconocimiento correcto de datos cuya etiqueta dentro del conjunto de datos es 1, es decir, hay abundancia de falsos negativos.
Ahora observemos una curva ROC para este modelo.
roc <- data.frame(
estimate = glm_test_syntetic_quadratic_class$estimate,
truth = glm_test_syntetic_quadratic_class$truth,
prob = predict_test_syntetic_quadratic)
yardstick::roc_curve(roc, truth, prob) %>%
autoplot() +
labs(title = "Curva ROC para modelo logístico")Ahora veamos el AUC para este modelo.
## # A tibble: 1 x 3
## .metric .estimator .estimate
## <chr> <chr> <dbl>
## 1 roc_auc binary 0.521
Note entonces que el ROC y el AUC son evidencia de que el modelo no tiene muy buena calidad, pues el AUC es cercano a 0.5 y además la curva ROC es muy cercana a la recta \(f(x) = x\), lo que demuestra que el clasificador es muy similar a un clasificador que etiqueta de manera aleatoria uniforme cada una de las filas del conjunto de datos.
e) Entrenamiento de un modelo con transformaciones no lineales
Ahora entrenaremos un modelo aplicando transformaciones no lineales sobre las características. En este caso aplicaremos las transformaciones \(X_1^2\) y \(X_2^2\).
nonlinear_glm_syntetic_quadratic <- glm(y ~ I(x1 ^ 2) + I(x2 ^ 2), data = syntetic_quadratic_train, family = "binomial")
summary(nonlinear_glm_syntetic_quadratic)##
## Call:
## glm(formula = y ~ I(x1^2) + I(x2^2), family = "binomial", data = syntetic_quadratic_train)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.403e-03 -2.000e-08 -2.000e-08 2.000e-08 1.590e-03
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -350.3 10118.8 -0.035 0.972
## I(x1^2) 23044.9 664686.0 0.035 0.972
## I(x2^2) -23048.1 665514.7 -0.035 0.972
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 5.1695e+02 on 374 degrees of freedom
## Residual deviance: 4.6856e-06 on 372 degrees of freedom
## AIC: 6
##
## Number of Fisher Scoring iterations: 25
Calculemos la tasa de clasificación correcta para este modelo.
predict_nonlinear_glm_syntetic_quadratic <- predict(nonlinear_glm_syntetic_quadratic, syntetic_quadratic_test, type = "response")
predict_class_nonlinear_glm_syntetic_quadratic <- rep(0, nrow(syntetic_quadratic_test))
predict_class_nonlinear_glm_syntetic_quadratic[predict_nonlinear_glm_syntetic_quadratic > 0.5] = 1
error_glm_nonlinear_syntetic_quadratic <- mean(predict_class_nonlinear_glm_syntetic_quadratic == syntetic_quadratic_test$y)
error_glm_nonlinear_syntetic_quadratic## [1] 1
Ahora, para realizar una análisis mas profundo, veamos la matriz de confusión para este modelo.
nonlinear_glm_syntetic_quadratic_class <- tibble(
pred = as.factor(predict_class_nonlinear_glm_syntetic_quadratic) ,
real = syntetic_quadratic_test$y)
nonlinear_glm_syntetic_quadratic_class %>% conf_mat(truth = real, estimate = pred) %>%
autoplot(type = "heatmap") +
labs(title = "Matriz de confusión para el modelo logístico (con transformaciones no lineales)")Observemos que la calidad de este modelo tiene un desempeño muy bueno, pues los residuales tienen una media de cero y por otro lado el AIC es menor, en comparación con los residuales de la regresión logística sin transformaciones cuadráticas. Esto demuestra entonces que cuando se tienen separaciones no lineales, se puede alcanzar una calidad muy buena en las predicciones sin recurrir a modelos sofisticados, y tomando en cuenta además que la regresión logística brinda mucha información adicional como por ejemplo, la importancia de las variables, el AIC y otros medidores que ayudan a comprender más a fondo el problema.
Miremos la curva roc para este modelo.
roc <- data.frame(
estimate = nonlinear_glm_syntetic_quadratic_class$pred,
truth = nonlinear_glm_syntetic_quadratic_class$real,
prob = predict_nonlinear_glm_syntetic_quadratic)
yardstick::roc_curve(roc, truth, prob) %>%
autoplot() +
labs(title = "Curva ROC para modelo logístico con transformaciones")Por otro lado, observemos el AUC para dicho modelo
## # A tibble: 1 x 3
## .metric .estimator .estimate
## <chr> <chr> <dbl>
## 1 roc_auc binary 1
Note que estas dos fuentes de información reafirman el hecho de que el modelo es muy bueno, pues la curva ROC esta lo mas alejada posible de la recta \(f(x) = x\), y su AUC es cercano a 1.
f) Gráficas sobre el modelo con transformaciones cuadráticas
Con el objetivo de reforzar el razonamiento presentado en el literal anterior, se presentará una gráfica de predicción sobre el conjunto de entrenamiento tomando en cuenta el modelo con las transformaciones sobre las variables propuestas anteriormente.
predict_nonlinear_train_syntetic_quadratic <- predict(nonlinear_glm_syntetic_quadratic, syntetic_quadratic_train, type = "response")
response_nonlinear_predict_train_syntetic_quadratic <- rep(0, nrow(syntetic_quadratic_train))
response_nonlinear_predict_train_syntetic_quadratic[predict_nonlinear_train_syntetic_quadratic > 0.5] = 1
ggplot(data = data.frame(x1 = syntetic_quadratic_train$x1,
x2 = syntetic_quadratic_train$x2,
y = as.factor(response_nonlinear_predict_train_syntetic_quadratic))) +
geom_point(mapping = aes(x = x1, y = x2, color = y), size = 2, alpha = 0.8) +
theme_economist() +
labs(x = TeX("X_1"), y = TeX("X_2"), color = "y", title = "Predicciones sobre el conjunto de entrenamiento (transformaciones cuadráticas)")g) Entrenamiento de un clasificador de soporte vectorial
Con el objetivo de realizar comparaciones con los métodos lineales generalizados presentados anteriormente, se entrenará un clasificador de soporte vectoríal sobre los datos sintéticos con separación no lineal, como se muestra a continuación.
svm_linear_syntetic_quadratic <- svm(y ~ ., data = syntetic_quadratic_train, kernel = "linear", cost = 10, scale = FALSE, probability = TRUE)
summary(svm_linear_syntetic_quadratic)##
## Call:
## svm(formula = y ~ ., data = syntetic_quadratic_train, kernel = "linear",
## cost = 10, probability = TRUE, scale = FALSE)
##
##
## Parameters:
## SVM-Type: C-classification
## SVM-Kernel: linear
## cost: 10
##
## Number of Support Vectors: 346
##
## ( 175 171 )
##
##
## Number of Classes: 2
##
## Levels:
## 0 1
Ahora realicemos una gráfica que muestre las fronteras de decisión para este modelo.
Note que el modelo no tiene la capacidad de evaluar correctamente la no linealidad dada para dichas clases. Como evidencia de ello, computaremos la tasa de clasificación correcta para este modelo sobre un conjunto de validación.
predict_test_svm_syntetic_quadratic <- predict(svm_linear_syntetic_quadratic, syntetic_quadratic_test, probability = TRUE)
error_linear_svm_syntetic_quadratic <- mean(predict_test_svm_syntetic_quadratic == syntetic_quadratic_test$y)
error_linear_svm_syntetic_quadratic## [1] 0.56
Ahora veamos cuál es la matriz de confusión para este modelo.
test_svm_syntetic_quadratic_class <- tibble(
pred = as.factor(predict_test_svm_syntetic_quadratic) ,
real = syntetic_quadratic_test$y)
test_svm_syntetic_quadratic_class %>% conf_mat(truth = real, estimate = pred) %>%
autoplot(type = "heatmap") +
labs(title = "Matriz de confusión para el clasificador de soporte vectorial")Note entonces que la calidad del modelo no es muy buena, debido a que es muy similar a tener un modelo de selección aleatoria. Es decir, si se tiene un modelo cuya respuesta es equiprobable para las clases 1 y 0 para toda entrada de \(X_1\) y \(X_2\), su desempeño va a ser muy similar al modelo presentado actualmente, lo cuál no es muy deseable. Esto se puede confirmar observando que la tasa de clasificación correcta de este modelo es 56%.
Para confirmar las afirmaciones presentadas en este punto, graficaremos una curva ROC con el fin de observar la calidad del modelo.
roc <- data.frame(
estimate = predict_test_svm_syntetic_quadratic,
truth = syntetic_quadratic_test$y,
prob = attr(predict_test_svm_syntetic_quadratic,"probabilities")[, "1"])
yardstick::roc_curve(roc, truth, prob) %>%
autoplot() +
labs(title = "Curva ROC para el SVC")Y extraemos su AUC de la siguiente manera.
## # A tibble: 1 x 3
## .metric .estimator .estimate
## <chr> <chr> <dbl>
## 1 roc_auc binary 0.511
Note que esto confirma lo dicho anteriormente, pues la curva ROC es muy cercana a la linea identidad y además su AUC es cercano a 0.5, lo que afirma que el clasificador es similar a un clasificador aleatorio
h) Entrenamiento de SVM con kernel no lineal
Ahora con el objetivo de mejorar los resultados sobre el clasificador de soporte vectorial, tomaremos una máquina de soporte vectorial con un kernel no lineal, en particular, tomaremos un kernel polinomial de grado dos, y lo compararemos con respecto al resultado del literal anterior.
svm_nonlinear_syntetic_quadratic <- svm(y ~ ., data = syntetic_quadratic_train, kernel = "polynomial", degree = 2, cost = 10, probability = TRUE)
summary(svm_nonlinear_syntetic_quadratic)##
## Call:
## svm(formula = y ~ ., data = syntetic_quadratic_train, kernel = "polynomial",
## degree = 2, cost = 10, probability = TRUE)
##
##
## Parameters:
## SVM-Type: C-classification
## SVM-Kernel: polynomial
## cost: 10
## degree: 2
## coef.0: 0
##
## Number of Support Vectors: 75
##
## ( 38 37 )
##
##
## Number of Classes: 2
##
## Levels:
## 0 1
Observemos una gráfica donde muestre las fronteras de decisión con el fin de evaluar la calidad del modelo.
Y ahora, calculemos la tasa de clasificación correcta para el modelo con kernel polinomial. Como se muestra a continuación
predict_svm_nonlinear_syntetic_quadratic <- predict(svm_nonlinear_syntetic_quadratic, syntetic_quadratic_test, probability = TRUE)
error_nonlinear_svm_syntetic_quadratic <- mean(predict_svm_nonlinear_syntetic_quadratic == syntetic_quadratic_test$y)
error_nonlinear_svm_syntetic_quadratic## [1] 0.96
Ahora si observamos la matriz de confusión de este modelo, se vera confirmada la calidad del modelo que se observó en la gráfica anterior
svm_nonlinear_syntetic_quadratic_class <- tibble(
pred = as.factor(predict_svm_nonlinear_syntetic_quadratic) ,
real = syntetic_quadratic_test$y)
svm_nonlinear_syntetic_quadratic_class %>% conf_mat(truth = real, estimate = pred) %>%
autoplot(type = "heatmap") +
labs(title = "Matriz de confusión para la máquina de soporte vectorial con kernel no lineal")Note entonces que la calidad de la predicción aumentó en gran medida, y no solo eso, además el modelo toma vectores de soporte mucho mas razonables que en el clasificador de soporte vectorial. Note además que la tasa de clasificación correcta es 96%, la cuál es una tasa de clasificación correcta extremadamente buena, obteniendo que el modelo representa correctamente el comportamiento de los datos.
Ahora para confirmar que la calidad del modelo es buena, grafiquemos la curva ROC.
roc <- data.frame(
estimate = predict_svm_nonlinear_syntetic_quadratic,
truth = syntetic_quadratic_test$y,
prob = attr(predict_svm_nonlinear_syntetic_quadratic,"probabilities")[, "1"])
yardstick::roc_curve(roc, truth, prob) %>%
autoplot() +
labs(title = "Curva ROC para la SVM con kernel no lineal")Note entonces que si computamos el AUC para este modelo se obtiene lo siguiente.
## # A tibble: 1 x 3
## .metric .estimator .estimate
## <chr> <chr> <dbl>
## 1 roc_auc binary 0.990
Se puede concluir a partir del AUC, que por una poca diferencia, el modelo logístico se desempeñó mucho mejor para este caso pues su valor es mas cercano a 1 que el valor mostrado anteriormente, aunque en la práctica no ofrece mucha diferencia ambas calificaciones.
i) Comentarios y conclusiones
Lo primero que se debe notar es que cuando los modelos no tienen una capacidad de representar una separación no lineal, como en el caso de la regresión logística tradicional y el clasificador de soporte vectorial, todo dependerá entonces del numero de elementos dentro de cada clase y su ubicación en el espacio. En general, cuando el número de elementos de cada clase es muy parejo y la distribución en el espacio es relativamente simétrica pero no linealmente separable, el clasificador obtiene un comportamiento de un clasificador aleatorio.
Por otro lado se observa que al dotar a un modelo de capacidad no lineal, este representa muy bien la frontera de decisión. Sin embargo es importante poner a juicio la pregunta: ¿en este caso es mejor un modelo logístico con transformaciones no lineales, o una máquina de soporte vectorial?, y en este caso se considera que es mas recomendable siempre abordar los problemas desde el modelo computacionalmente más sencillo, el cuál es, el modelo logístico; además es importante resaltar que el modelo logístico brinda mucha más información respecto a los predictores que la máquina de soporte vectorial, sin dejar de lado la simplicidad e interpretabilidad del modelo, pues es importante recordar que este modelo ofrece una interpretabilidad probabilística, y da la posibilidad de abordar los resultados desde un punto de vista mucho mas estadístico en comparación con el SVM.
Sección 9.7 - Ejercicio 6
a) Generar un conjunto de datos apenas linealmente separable
Se crea un conjunto de datos apenas linealmente separable Con el propósito de probar que una máquina de soporte vectorial con un valor de costo que clasifica de manera incorrecta en el conjunto de entrenamiento utilizando un valor de costo pequeño puede llegar a ser mejor en el conjunto de prueba que uno con costo mejor no cometa ningún error en el conjunto de entrenamiento.
Entendemos por conjunto apenas linealmente separable un conjunto tal que sobre la línea de separación que utilizará el clasificador, haya observaciones de ambas clases.
set.seed(1)
synth <- tibble(a = c(runif(200, min = 0, max = 2), runif(200, min = 0, max = 2)),
b = c(runif(200, min = 0, max = 1.1), runif(200, min = 1, max = 2)),
resp = as_factor(c(rep("clase 1", 200), rep("clase 2", 200))))
synth %>%
ggplot(aes(x = a, y = b, color = resp)) +
geom_point()## # A tibble: 6 x 3
## a b resp
## <dbl> <dbl> <fct>
## 1 0.531 0.725 clase 1
## 2 0.744 0.204 clase 1
## 3 1.15 1.05 clase 1
## 4 1.82 0.988 clase 1
## 5 0.403 1.04 clase 1
## 6 1.80 0.796 clase 1
b) Cómputo de los errores de validación utilizando SVM con diferentes rangos de costo
tune_synth <- tune(svm, resp ~ ., data = synth, kernel = "linear",
ranges = list(cost = seq(from = 0.01, to = 50, by = 0.1)))
l <- ggplot(data = tune_synth$performances) +
geom_line(mapping = aes(x = cost, y = error)) +
scale_x_log10() +
theme_economist() +
labs(title = "Errores respecto al costo", x = "Costo", y = "Error")
ggplotly(l)##
## Call:
## best.tune(method = svm, train.x = resp ~ ., data = synth, ranges = list(cost = seq(from = 0.01,
## to = 50, by = 0.1)), kernel = "linear")
##
##
## Parameters:
## SVM-Type: C-classification
## SVM-Kernel: linear
## cost: 9.01
##
## Number of Support Vectors: 48
##
## ( 24 24 )
##
##
## Number of Classes: 2
##
## Levels:
## clase 1 clase 2
Se identifica que utilizando validación cruzada, se cometen errores considerablemente bajos (0.04) con un costo entre 0.21 y 1, sin embargo, los errores de clasificación incrementan cuando el costo está entre 2.71 y 50;lo cual indicaría que no necesariamente aumentar el costo a un valor alto mejora la predicción en el conjunto de prueba.
c) Creación de un conjunto de validación y valores de costo asociados
split_ratio_synth <- 0.75
smp_size_synth <- floor(split_ratio_synth * nrow(synth))
train_ind_synth <- sample(seq_len(nrow(synth)), size = smp_size_synth)
synth_train <- synth[train_ind_synth, ]
synth_test <- synth[-train_ind_synth, ]
dim(synth_train)## [1] 300 3
## [1] 100 3
Modelo SVM con costo seleccionado por validación cruzada
Gráfica de clasificación en el conjunto de entrenamiento:
##
## Call:
## best.tune(method = svm, train.x = resp ~ ., data = synth, ranges = list(cost = seq(from = 0.01,
## to = 50, by = 0.1)), kernel = "linear")
##
##
## Parameters:
## SVM-Type: C-classification
## SVM-Kernel: linear
## cost: 9.01
##
## Number of Support Vectors: 48
##
## ( 24 24 )
##
##
## Number of Classes: 2
##
## Levels:
## clase 1 clase 2
Clasificación en el conjunto de prueba:
En general se evidencia que la máquina de soporte vectorial con kernel lineal no se vio muy afectada al no conocer la totalidad de los datos del conjunto.
svm_synth_class <- tibble(
pred = as.factor(predict(tune_synth$best.model, synth_test)) ,
real = synth_test$resp)
svm_synth_class %>% conf_mat(truth = real, estimate = pred) %>%
autoplot(type = "heatmap") +
labs(title = "Matriz de confusión para la SVM con costo óptimo por validación cruzada")En este caso, la máquina de soporte vectorial con kernel lineal tiene una clasificación casi perfecta en el conjunto de prueba, con una predicción incorrecta de la clase 1 cuando en realidad era la clase 2, ahora validemos la máquina de soporte vectorial con costo elevado.
Modelo con costo = 1000
svm_synth_50 <- svm(resp ~ ., data = synth_train, kernel = "linear", cost = 1000, scale = FALSE, probability = TRUE)
summary(svm_synth_50)##
## Call:
## svm(formula = resp ~ ., data = synth_train, kernel = "linear", cost = 1000,
## probability = TRUE, scale = FALSE)
##
##
## Parameters:
## SVM-Type: C-classification
## SVM-Kernel: linear
## cost: 1000
##
## Number of Support Vectors: 29
##
## ( 14 15 )
##
##
## Number of Classes: 2
##
## Levels:
## clase 1 clase 2
Gráfica de clasificación en el conjunto de entrenamiento
Clasificación en el conjunto de prueba:
svm_synth_class_50 <- tibble(
pred = as.factor(predict(svm_synth_50, synth_test, probability = TRUE)) ,
real = synth_test$resp)
svm_synth_class_50 %>% conf_mat(truth = real, estimate = pred) %>%
autoplot(type = "heatmap") +
labs(title = "Matriz de confusión para la SVM con costo = 1000")En este caso, la máquina de soporte vectorial con kernel lineal y costo alto, tiene también una clasificación buena en el conjunto de prueba, sin embargo esta es menor que en el modelo con costo óptimo. Esto es un buen indicador de que aumentar el costo no necesariamente incrementa la capacidad predictiva de un modelo en el conjunto de prueba.
Sección 9.7 - Ejercicio 7
En este caso, usaremos un SVM con el fin de predecir en qué momento un carro tiene alto o bajo consumo de combustible en el conjunto de datos Auto, considerando cuándo el consumo está por encima o por debajo de la mediana.
a) Preparación del conjunto de datos
Inicialmente crearemos una variable binaria que tomará el valor de 1 cuando el consumo de combustible está por encima de la media y tomará un valor de 0 en caso contrario.
auto <- as_tibble(Auto)
median_mpg <- median(auto$mpg)
auto <- auto %>% mutate(
mpg01 = case_when(
mpg >= median_mpg ~ 1,
mpg < median_mpg ~ 0
),
mpg01 = as.factor(mpg01),
origin = as.factor(origin)
)
# Retiramos la variable mpg
auto_fit <- auto %>% dplyr::select(-mpg, -name)
head(auto)## # A tibble: 6 x 10
## mpg cylinders displacement horsepower weight acceleration year origin name
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <fct> <fct>
## 1 18 8 307 130 3504 12 70 1 chev~
## 2 15 8 350 165 3693 11.5 70 1 buic~
## 3 18 8 318 150 3436 11 70 1 plym~
## 4 16 8 304 150 3433 12 70 1 amc ~
## 5 17 8 302 140 3449 10.5 70 1 ford~
## 6 15 8 429 198 4341 10 70 1 ford~
## # ... with 1 more variable: mpg01 <fct>
Ahora observemos una gráfica descriptiva de los datos.
Al notar la imágen anterior, se puede observar que muchas de las variables tienen una separación notoria. Si las separaciones son lineales, es posible trazar hiperplanos que separen los datos de manera adecuada, sobre todo en las variables displacement, weight y cilinders
b) Clasificador de soporte vectorial
Construyamos un clasificador de soporte vectorial y observemos varios valores para la varible cost usando validación cruzada, para ello consideremos el siguiente codigo, y con el se extraerá el mejor clasificador con respecto a la variable cost.
set.seed(20)
tune_svc <- tune(svm, mpg01 ~ ., data = auto_fit, kernel = "linear",
ranges = list(cost = seq(from = 0.01, to = 5, by = 0.1)))
summary(tune_svc)##
## Parameter tuning of 'svm':
##
## - sampling method: 10-fold cross validation
##
## - best parameters:
## cost
## 1.41
##
## - best performance: 0.08429487
##
## - Detailed performance results:
## cost error dispersion
## 1 0.01 0.08685897 0.04728478
## 2 0.11 0.09955128 0.04897224
## 3 0.21 0.09698718 0.04946574
## 4 0.31 0.09692308 0.03768179
## 5 0.41 0.10205128 0.04515964
## 6 0.51 0.09185897 0.03245902
## 7 0.61 0.08685897 0.03266541
## 8 0.71 0.08685897 0.03266541
## 9 0.81 0.08942308 0.03699143
## 10 0.91 0.08685897 0.03483004
## 11 1.01 0.08685897 0.03483004
## 12 1.11 0.08685897 0.03483004
## 13 1.21 0.08685897 0.03483004
## 14 1.31 0.08685897 0.03483004
## 15 1.41 0.08429487 0.03448754
## 16 1.51 0.08429487 0.03448754
## 17 1.61 0.08429487 0.03448754
## 18 1.71 0.08429487 0.03448754
## 19 1.81 0.08429487 0.03448754
## 20 1.91 0.08429487 0.03448754
## 21 2.01 0.08429487 0.03448754
## 22 2.11 0.08429487 0.03448754
## 23 2.21 0.08429487 0.03448754
## 24 2.31 0.08429487 0.03448754
## 25 2.41 0.08429487 0.03448754
## 26 2.51 0.08429487 0.03448754
## 27 2.61 0.08429487 0.03448754
## 28 2.71 0.08429487 0.03448754
## 29 2.81 0.08429487 0.03448754
## 30 2.91 0.08429487 0.03448754
## 31 3.01 0.08429487 0.03448754
## 32 3.11 0.08429487 0.03448754
## 33 3.21 0.08429487 0.03448754
## 34 3.31 0.08685897 0.03686780
## 35 3.41 0.08685897 0.03686780
## 36 3.51 0.08685897 0.03686780
## 37 3.61 0.08685897 0.03686780
## 38 3.71 0.08685897 0.03686780
## 39 3.81 0.08685897 0.03686780
## 40 3.91 0.08685897 0.03686780
## 41 4.01 0.08685897 0.03686780
## 42 4.11 0.08685897 0.03686780
## 43 4.21 0.08685897 0.03686780
## 44 4.31 0.08685897 0.03686780
## 45 4.41 0.08685897 0.03686780
## 46 4.51 0.08685897 0.03686780
## 47 4.61 0.08685897 0.03686780
## 48 4.71 0.08685897 0.03686780
## 49 4.81 0.08685897 0.03686780
## 50 4.91 0.08685897 0.03686780
Notemos entonces que el mejor modelo corresponde a aquel que posee un costo de 1.41, obteniendo asi un error de 0.0842949. Sin embargo, es interesante observar como se comporta el error respecto a la variable costo de manera gráfica para los valores probados.
g <- ggplot(data = tune_svc$performances) +
geom_line(mapping = aes(x = cost, y = error)) +
theme_economist() +
labs(title = "Errores respecto al costo", x = "Costo", y = "Error")
ggplotly(g)Note pues que se pudo tomar cualquier costo entre 1.41 y 3.21 y se hubiera obtenido el mismo resultado.
El desempeño del modelo tiene una buena calidad, pues reporta una tasa de clasificación incorrecta de 8.4294872%. Es posible que este resultado haya sido debido a que en la gráfica descriptiva, se noto una clara separación de múltiples variables al resaltarlas para cada una de las clases, y es posible que el modelo haya detectado y expresado de manera correcta dicha separación.
c) SVM con kernel no lineal
Ahora encontremos un modelo óptimo para este conjunto de datos usando un modelo SVM con kernel no lineal. Primero, usaremos un modelo con un kernel polinomial y extraeremos el mejor modelo probando sobre varios grados del polinomio; posteriormente, probaremos con una máquina de soporte vectorial con kernel radial y encontraremos el mejor modelo para varios valores de \(\gamma\). En ambos modelos, se buscará un valor adecuado para el costo.
SVM con kernel polinomial
A continuación entrenaremos un modelo con kernel polinomial y buscaremos el mejor modelo usando validación cruzada para encontrar el mejor valor para el grado del polinomio.
set.seed(20)
tune_svm_polynomial <- tune(svm, mpg01 ~ ., data = auto_fit, kernel = "polynomial", probability = TRUE,
ranges = list(cost = seq(from = 0.01, to = 5, by = 0.2),
degree = seq(from = 2, to = 10, by = 1)))Se obtiene que el mejor modelo tiene un grado 3 y un costo de 3.01, para así obtener un tasa de error de clasificación de 7.4166667%, el cual es levemente mejor que la máquina de soporte vectorial con kernel lineal.
Realicemos una gráfica para observar el comportamiento del error respecto al costo y al grado
fig <- plot_ly(data = tune_svm_polynomial$performances, x=~cost,y=~degree, z=~ -1/(error ^ 2), type = "contour", colorscale='Jet') %>%
layout(
xaxis = list(title = "Costo"),
yaxis = list(title = "Grado del polinomio"),
title = "Gráfica sobre el ")
fig %>% colorbar(title = "-1/(Error ^ 2)")Nota: En esta figura no se graficó el error \(E\) neto como lo presentan los modelos, pues en este caso la gráfica no expresa de manera clara los diferentes niveles para algunas regiones, por tanto se gráfica en realidad \(-1/E^2\), con el fin de observar de mejor manera el comportamiento del error en este caso.
Esta figura se puede observar que es una mala estrategia tomar valores para el grado del polinomio mucho más altos, aunque se pensaría que estos expresan una variedad mayor de funciones porque tienen muchos grados de libertad, es precisamente esto lo que causa que el modelo se sobreajuste a los datos de entrenamiento, por tanto, en este caso es mejor tomar valores del grado del polinomio mucho mas bajos para obtener un modelo que logre ser lo más general posible para los datos considerados.
SVM con kernel radial
Ahora, entrenemos una máquina de soporte vectorial con un kernel radial. Para ello encontraremos los mejores parámetros para el costo y para \(\gamma\).
set.seed(20)
tune_svm_radial <- tune(svm, mpg01 ~ ., data = auto_fit, kernel = "radial",
ranges = list(cost = seq(from = 0.01, to = 5, by = 0.3),
gamma = seq(from = 0.01, to = 1, by = 0.02)))En este caso, los mejores parámetros encontrados que minimizan la tasa de error de clasificación corresponden a \(\gamma\) con un valor de 0.53 y un costo de 1.21, para así obtener una tasa de clasificación erronea de 6.1346154%.
Ahora, con el objetivo de visualizar como varía el error de acuerdo a \(\gamma\) y el costo, veamos la siguiente gráfica.
fig <- plot_ly(data = tune_svm_radial$performances, x=~cost,y=~gamma, z=~ -1/(error ^ 2), type = "contour", colorscale='Jet') %>%
layout(
xaxis = list(title = "Costo"),
yaxis = list(title = "Gamma"))
fig %>% colorbar(title = "-1/(Error ^ 2)")Nota: En esta figura no se graficó el error \(E\) neto como lo presentan los modelos, pues en este caso la gráfica no expresa de manera clara los diferentes niveles para algunas regiones, por tanto se gráfica en realidad \(-1/E^2\), con el fin de observar de mejor manera el comportamiento del error en este caso.
Esta gráfica demuestra que el modelo es altamente sensible al costo, al contrario del parrametro \(\gamma\) que no ofrece una variabilidad tan alta. Es decir, note que el error solo cambia en gran medida si el costo es demasiado pequeño; si se toma un costo muy pequeño, se obtendrá un error alto sin importar si \(\gamma\) se mueve en un rango bastante amplio. Otra conclusión que se puede extraer es que si tomamos un costo mayor a 0.31, y seleccionamos algún valor para \(\gamma\) entre 0 y 1, se obtienen valores del error con un cambio pequeño.
d) Gráficas
De acuerdo a la gráfica de correlaciones presentada al inicio de la sección, se puede observar que la interacción de los caballos de fuerza del automóvil tienen una interacción notoria con las demás variables y a su vez muestra una tendencia marcada al diferenciar los datos para cada una de las clases de la variable mpg01. En vista de ello, centraremos nuestro análisis en tomar la variable horsepowery compararla con las variables que ofrecen más información.
Para este análisis comprararemos la variable horsepower con otras tres variables, las cuáles son acceleration, displacement y weight, observando allí como actua cada uno de los modelos. Es importante aclarar que en cada una de las gráficas se tomo una tajada (o slice) justo en la mediana de cada una de las variables restantes.
SVM con kernel lineal
Inicialmente si observemos las gráficas de las máquinas de soporte vectorial para cada pareja de variables. A continuación mostraremos la gráfica de la máquina de soporte vectorial para las variables horsepower y acceleration.
plot(tune_svc$best.model, auto_fit, horsepower ~ acceleration,
slice = list(year = median(auto_fit$year),
cylinders = median(auto_fit$cylinders),
displacement = median(auto_fit$displacement),
weight = median(auto_fit$weight),
origin = 1))En este caso, note que la máquina de soporte vectorial reconoce muy bien la separación mostrada en este par de variables, por tanto la máquina de soporte vectorial con kernel lineal se muestra efectiva en esta situación.
Ahora grafiquemos las variables horsepower y displacement y observemos la frontera de decisión en este caso.
plot(tune_svc$best.model, auto_fit, horsepower ~ displacement,
slice = list(year = median(auto_fit$year),
cylinders = median(auto_fit$cylinders),
acceleration = median(auto_fit$acceleration),
weight = median(auto_fit$weight),
origin = 1))Se encuentra entonces que hay ciertos carros cuyo consumo de combustible está por debajo de la mediana pero aún así aparecen clasificados como autos con alto consumo. Esto puede ocurrir debido al gran peso estadístico que ofrecen los datos para la clase de los automóviles de alto consumo.
Finalmente para este modelo, grafiquemos la frontera de decisión para las variables horsepower y weight.
plot(tune_svc$best.model, auto_fit, horsepower ~ weight,
slice = list(year = median(auto_fit$year),
acceleration = median(auto_fit$acceleration),
displacement = median(auto_fit$displacement),
cylinders = median(auto_fit$cylinders),
origin = 1))Este caso es interesante, pues note que la separación lineal es razonablemente adecuada considerando que estamos en un caso no separable, y aún así el modelo se comporta de la manera esperada clasificando correctamente a la mayoría de los datos.
SVM con kernel polinomial
En esta sección, consideraremos las mismas parejas que se mostraron en el análisis anterior para una máquina de soporte vectorial con un kernel polinomial de grado 3.
Para las variables horsepower y acceleration se tiene la siguiente frontera de decisión.
plot(tune_svm_polynomial$best.model, auto_fit, horsepower ~ acceleration,
slice = list(year = median(auto_fit$year),
cylinders = median(auto_fit$cylinders),
displacement = median(auto_fit$displacement),
weight = median(auto_fit$weight),
origin = 1))Note que en este caso, la clase de los automóviles con mayor consumo de combustible es la menos favorecida, pues gran parte de automóviles en esta clase fue considerada como automóviles de bajo consumo. Sin embargo la frontera de decisión representa adecuadamente la forma en la que se distribuyen los autos de mayor consumo. Es importante resaltar también que esta frontera de decisión, a pesar de ser basada en un kernel polinomial de grado mayor a uno, parece como si fuera linel.
Ahora grafiquemos las frontera de decisión para las variables horsepower y weight.
plot(tune_svm_polynomial$best.model, auto_fit, horsepower ~ weight,
slice = list(year = median(auto_fit$year),
acceleration = median(auto_fit$acceleration),
displacement = median(auto_fit$displacement),
cylinders = median(auto_fit$cylinders),
origin = 1))Una vez más, la frontera separa de una manera adecuada las clases en esta situación. Sin embargo, note que la frontera de decisión parece lineal, así como en el caso anterior.
Finalmente grafiquemos la frontera de decisión para la variable horsepower y displacement.
plot(tune_svm_polynomial$best.model, auto_fit, horsepower ~ displacement,
slice = list(year = median(auto_fit$year),
cylinders = median(auto_fit$cylinders),
acceleration = median(auto_fit$acceleration),
weight = median(auto_fit$weight),
origin = 1))Note que en este caso la separación es efectiva, sin embargo hay algunos datos que no permiten que el conjunto sea separable, lo que no da oportunidad al clasificador de diferenciar estos casos particulares. Aun así, el clasificador etiqueta de manera adecuada los datos de entrenamiento segun su tendencia en la mayoría de los casos. Como sucedió en los demás casos, esta frontera no tiene una diferencia marcada en comparación a una frontera producida por un kernel lineal.
SVM con kernel radial
En este caso, graficaremos las variables anteriormente mencionadas sobre un clasificador con kernel radial y observaremos el desempeño del calsificador en algunas de las situaciones.
Inicialmente veamos la gráfica para las variables horsepower y acceleration.
plot(tune_svm_radial$best.model, auto_fit, horsepower ~ acceleration,
slice = list(year = median(auto_fit$year),
cylinders = median(auto_fit$cylinders),
displacement = median(auto_fit$displacement),
weight = median(auto_fit$weight),
origin = 1))Note que para esta sección, se muestra claramente una frontera no lineal, sin embargo, a este punto no hay una separación notablemente efectiva. Esto pudo ocurrir probablemente al tomar una sección que no es adecuada para las demas variables. En esta gráfica, muchos de los autos con bajo consumo son clasificados como autos de alto consumo, pero esto puede corregirse si se toma una tajada más adecuada.
Con la siguiente gráfica vamos a visualizar las variables horsepower y weight.
plot(tune_svm_radial$best.model, auto_fit, horsepower ~ weight,
slice = list(year = median(auto_fit$year),
acceleration = median(auto_fit$acceleration),
displacement = median(auto_fit$displacement),
cylinders = median(auto_fit$cylinders),
origin = 1))Note que en este caso se nota también un el uso de un kernel no lineal por la forma de la frontera de decisión a diferencia del kernel polinomial. Se nota entonces que la separación le da mucha prioridad a los casos positivos, es decir, en esta tajada, la frontera trata de abarcar todos los carros con alto consumo, pues probablemente en esta sección haya una gran cantidad de datos ubicados en esta tajada que estan clasificados como 1.
Finalmente, grafiquemos la frontera de decisión para las variables horsepower y displacement.
plot(tune_svm_radial$best.model, auto_fit, horsepower ~ displacement,
slice = list(year = median(auto_fit$year),
cylinders = median(auto_fit$cylinders),
acceleration = median(auto_fit$acceleration),
weight = median(auto_fit$weight),
origin = 1))En este caso ocurre lo mismo que el anterior caso, la frontera de decisión trata de abarcar la mayor cantidad de automóviles de alto consumo. Cabe aclarar que en estos casos no quiere decir que la calidad de la clasificación sea mala, por el contrario, puede suceder que en esta tajada hayan muchos casos en los cuales los automóviles tengan alto consumo y la frontera de decisión tiene que ampliarse para acaparar a la mayor cantidad de casos posibles.
Curvas ROC y AUC
Con el fin de comparar los modelos de clasificación, vamos a graficar las curvas ROC y a calcular el AUC para cada uno de los modelos usando un conjunto de entrenamiento y otro de validación, y así observar qué modelo se desempeña mejor en este caso.
Inicialmente realicemos una partición de los datos en un conjunto de entrenamiento y otro de validación.
split_ratio <- 0.75
smp_size <- floor(split_ratio * nrow(auto_fit))
train_ind <- sample(seq_len(nrow(auto_fit)), size = smp_size)
auto_fit_train <- auto_fit[train_ind, ]
auto_fit_test <- auto_fit[-train_ind, ]Entrenemos un clasificador de soporte vectorial con los parámetros óptimos obtenidos en literales anteriores y tabulemos los resultados.
svm_linear <- svm(mpg01 ~ .,
data = auto_fit_train,
cost = tune_svc$best.parameters$cost,
kernel = "linear",
probability = TRUE)
svm_linear_test <- predict(svm_linear, auto_fit_test, probability = TRUE)
roc_linear <- data.frame(
estimate = svm_linear_test,
truth = auto_fit_test$mpg01,
prob = attr(svm_linear_test,"probabilities")[, "1"],
Modelo = "Lineal")
head(roc_linear)## estimate truth prob Modelo
## 1 0 0 1.716515e-03 Lineal
## 2 0 0 9.973473e-04 Lineal
## 3 0 0 7.900009e-04 Lineal
## 4 1 1 9.599086e-01 Lineal
## 5 0 0 2.763476e-06 Lineal
## 6 0 0 8.960600e-06 Lineal
Ahora, entrenemos una máquina de soporte vectorial con kernel polinomial, usando el valor del grado óptimo encontrado en los literales anteriores.
svm_poly <- svm(mpg01 ~ .,
data = auto_fit_train,
cost = tune_svm_polynomial$best.parameters$cost,
degree = tune_svm_polynomial$best.parameters$degree,
kernel = "polynomial",
probability = TRUE)
svm_poly_test <- predict(svm_poly, auto_fit_test, probability = TRUE)
roc_poly <- data.frame(
estimate = svm_poly_test,
truth = auto_fit_test$mpg01,
prob = attr(svm_poly_test,"probabilities")[, "1"],
Modelo = "Polinomial")
head(roc_poly)## estimate truth prob Modelo
## 1 0 0 0.0000103496 Polinomial
## 2 0 0 0.0000001000 Polinomial
## 3 0 0 0.0000001000 Polinomial
## 4 1 1 0.9972020683 Polinomial
## 5 0 0 0.0000001000 Polinomial
## 6 0 0 0.0000001000 Polinomial
Finalmente entrenemos una máquina de soporte vectorial con kernel radial usando el valor de \(\gamma\) encontrado en los literales anteriores.
svm_radial <- svm(mpg01 ~ .,
data = auto_fit_train,
cost = tune_svm_radial$best.parameters$cost,
gamma = tune_svm_radial$best.parameters$gamma,
kernel = "radial",
probability = TRUE)
svm_radial_test <- predict(svm_radial, auto_fit_test, probability = TRUE)
roc_radial <- data.frame(
estimate = svm_radial_test,
truth = auto_fit_test$mpg01,
prob = attr(svm_radial_test,"probabilities")[, "1"],
Modelo = "Radial")
head(roc_radial)## estimate truth prob Modelo
## 1 0 0 0.04368034 Radial
## 2 0 0 0.11574914 Radial
## 3 0 0 0.08710994 Radial
## 4 1 1 0.92136184 Radial
## 5 0 0 0.05301847 Radial
## 6 0 0 0.04964640 Radial
Teniendo la información recolectada anteriormente, grafiquemos las curvas ROC para cada uno de los modelos usados en este caso.
roc_merge <- rbind(roc_linear, roc_poly, roc_radial)
roc_merge %>%
group_by(Modelo) %>%
roc_curve(truth, prob) %>%
autoplot() +
labs(title = "Curvas ROC")Notemos entonces que a partir de las curvas ROC no se puede afirmar que hay un modelo predominante. Se puede notar que los modelos entrenados tienen una buena calidad pues las curvas ROC estan alejadas de la recta dada por la función \(f(x) = x\).
Con el objetivo de determinar una diferencia entre los modelos, veamos el AUC para cada uno de ellos.
## # A tibble: 3 x 4
## Modelo .metric .estimator .estimate
## <fct> <chr> <chr> <dbl>
## 1 Lineal roc_auc binary 0.972
## 2 Polinomial roc_auc binary 0.966
## 3 Radial roc_auc binary 0.956
En este caso, el modelo lineal tiene un mejor desempeño, pues el área bajo la curva ROC (AUC) para dicho modelo es mayor a los demas modelos. Sin embargo es importante resaltar que los modelos no tienen una diferencia muy marcada entre ellos pues el valor del AUC no varía demasiado, aun así, el modelo con kernel lineal posee mayor simplicidad y obtiene un mayor AUC, por tanto es la mejor opción para realizar predicciones sobre este conjunto de datos.
Sección 9.7 - Ejercicio 8
a) Conjuntos de entrenamiento y validación
El conjunto de datos OJ contiene 1070 compras donde los clientes compraron un jugo de naranja de Citrus Hill o de Maid Orange Juice y se guardaron un número de características de los clientes. Se crea un conjunto de entrenamiento con 800 observaciones y 270 de prueba.
oj_train <- as_tibble(OJ) %>%
rownames_to_column() %>%
sample_n(size = 800)
oj_test <- as_tibble(OJ) %>%
rownames_to_column() %>%
anti_join(oj_train, by = "rowname")
oj_train <- oj_train %>% dplyr::select(-rowname)
oj_test <- oj_test %>% dplyr::select(-rowname)
dim(oj_train)## [1] 800 18
## [1] 270 18
b) Máquina de soporte vectorial con costo = 0.01 para clasificar el tipo de jugo comprado
Dado que la cantidad de variables del conjunto OJ es relativamente grande, la flexibilidad o complejidad añadida por un kernel no lineal debería ser innecesaria, por este motivo se utilizará un kernel lineal en este apartado y en los apartados f) y g) verificaremos los resultados con kernel radial y polinomial.
svm_oj <- svm(Purchase ~ ., data = oj_train, kernel = "linear", cost = 0.01, probability = TRUE)
summary(svm_oj)##
## Call:
## svm(formula = Purchase ~ ., data = oj_train, kernel = "linear", cost = 0.01,
## probability = TRUE)
##
##
## Parameters:
## SVM-Type: C-classification
## SVM-Kernel: linear
## cost: 0.01
##
## Number of Support Vectors: 445
##
## ( 224 221 )
##
##
## Number of Classes: 2
##
## Levels:
## CH MM
En total se tienen 441 vectores de soporte, de los cuales 220 corresponden al jugo Citru Hill y 221 al jugo Maid Orange Juice.
c) Errores de entrenamiento y validación
training_oj_class <- tibble(pred = svm_oj$fitted,
real = oj_train$Purchase)
training_oj_class %>% conf_mat(truth = real, estimate = pred) %>%
autoplot(type = "heatmap") +
labs(title = "Matriz de confusión de la SVM con costo = 0.01 en el conjunto de entrenamiento")## # A tibble: 2 x 3
## .metric .estimator .estimate
## <chr> <chr> <dbl>
## 1 accuracy binary 0.828
## 2 kap binary 0.640
En este caso, se tiene una tasa de clasificación correcta del 82.2% en el conjunto de entrenamiento, pero esta métrica no es de tanto interés como la del conjunto de validación, ya que se quiere probar el rendimiento de este algoritmo con información desconocida para el modelo.
testing_oj_class <- tibble(pred = predict(svm_oj, newdata = oj_test),
real = oj_test$Purchase)
testing_oj_class %>% conf_mat(truth = real, estimate = pred) %>%
autoplot(type = "heatmap") +
labs(title = "Matriz de confusión de la SVM con costo = 0.01 en el conjunto de prueba")## # A tibble: 2 x 3
## .metric .estimator .estimate
## <chr> <chr> <dbl>
## 1 accuracy binary 0.852
## 2 kap binary 0.675
De manera atípica, el algoritmo tiene mejor rendimiento en el conjunto de prueba, con una tasa de clasificación correcta del 83.3%.
d) Validación cruzada para hacer tunning de la mejor SVM
tune_oj <- tune(svm, Purchase ~ ., data = oj_train, kernel = "linear", probability = TRUE,
ranges = list(cost = seq(from = 0.01, to = 10, by = 1)))
k <- ggplot(data = tune_oj$performances) +
geom_line(mapping = aes(x = cost, y = error)) +
scale_x_log10() +
theme_economist() +
labs(title = "Errores respecto al costo", x = "Costo", y = "Error")
ggplotly(k)##
## Call:
## best.tune(method = svm, train.x = Purchase ~ ., data = oj_train,
## ranges = list(cost = seq(from = 0.01, to = 10, by = 1)), kernel = "linear",
## probability = TRUE)
##
##
## Parameters:
## SVM-Type: C-classification
## SVM-Kernel: linear
## cost: 9.01
##
## Number of Support Vectors: 335
##
## ( 167 168 )
##
##
## Number of Classes: 2
##
## Levels:
## CH MM
El mejor modelo utiliza 171 vectores de soporte para la clase Citru Hill y 172 para el jugo Maid Orange Juice utilizando un costo de 2.16, considerablemente menos vectores que el modelo con costo = 0.01 que tenía un error de validación cruzada de 0.17 como se puede ver en la gráfica.
e) Errores de entrenamiento y validación para el mejor modelo seleccionado
training_oj_lineal_class <- tibble(pred = tune_oj$best.model$fitted,
real = oj_train$Purchase)
training_oj_lineal_class %>% conf_mat(truth = real, estimate = pred) %>%
autoplot(type = "heatmap") +
labs(title = "Matriz de confusión de la SVM con costo = 2.16 en el conjunto de entrenamiento")## # A tibble: 2 x 3
## .metric .estimator .estimate
## <chr> <chr> <dbl>
## 1 accuracy binary 0.842
## 2 kap binary 0.672
Se tiene una tasa de clasificación correcta del 83.1% en el conjunto de entrenamiento, una mejora ligera con respecto al clasificador con costo = 0.01, sin embargo, esta mejora es en entrenamiento y no representa el comportamiento del algoritmo para predecir de manera correcta. Realicemos la validación en el conjunto de prueba.
testing_oj_lineal_class <- tibble(pred = predict(tune_oj$best.model, newdata = oj_test),
real = oj_test$Purchase)
testing_oj_lineal_class %>% conf_mat(truth = real, estimate = pred) %>%
autoplot(type = "heatmap") +
labs(title = "Matriz de confusión de la SVM con costo = 0.01 en el conjunto de prueba")## # A tibble: 2 x 3
## .metric .estimator .estimate
## <chr> <chr> <dbl>
## 1 accuracy binary 0.848
## 2 kap binary 0.668
La mejora utilizando el mejor costo por validación cruzada es muy ligera, cambiando de 83.3% a 83.7%, lo cual puede ser significativo para la decisión que se toma con respecto a una cantidad grande de jugos de naranja.
d y e) Comparación de los modelos con kernel radial y polinómico de grado 2
Ajustemos inicialmente los 2 nuevos modelos con el conjunto de entrenamiento utilizando la función tune para seleccionar el mejor modelo directamente.
tune_oj_radial <- tune(svm, Purchase ~ ., data = oj_train, kernel = "radial", probability = TRUE,
ranges = list(cost = seq(from = 0.01, to = 10, by = 0.1)))
tune_oj_poly <- tune(svm, Purchase ~ ., data = oj_train, kernel = "polynomial", degree = 2, probability = TRUE,
ranges = list(cost = seq(from = 0.01, to = 10, by = 0.1)))
summary(tune_oj_radial$best.model)##
## Call:
## best.tune(method = svm, train.x = Purchase ~ ., data = oj_train,
## ranges = list(cost = seq(from = 0.01, to = 10, by = 0.1)), kernel = "radial",
## probability = TRUE)
##
##
## Parameters:
## SVM-Type: C-classification
## SVM-Kernel: radial
## cost: 0.91
##
## Number of Support Vectors: 380
##
## ( 187 193 )
##
##
## Number of Classes: 2
##
## Levels:
## CH MM
##
## Call:
## best.tune(method = svm, train.x = Purchase ~ ., data = oj_train,
## ranges = list(cost = seq(from = 0.01, to = 10, by = 0.1)), kernel = "polynomial",
## degree = 2, probability = TRUE)
##
##
## Parameters:
## SVM-Type: C-classification
## SVM-Kernel: polynomial
## cost: 3.81
## degree: 2
## coef.0: 0
##
## Number of Support Vectors: 380
##
## ( 188 192 )
##
##
## Number of Classes: 2
##
## Levels:
## CH MM
Para el modelo con kernel polinomial se tienen 372 vectores de soporte y un costo de 5.01; para el modelo con kernel radial se tienen 407 vectores de soporte y un costo de 0.61.
Comparación de los 3 modelos de SVM en el conjunto de prueba
Para comparar efectivamente estos modelos sin tener que analizar 6 matrices de confusión, utilizaremos la tasa de verdaderos positivos (sensibilidad) vs la tasa de falsos positivos (1 - especificidad) mediante una curva ROC en el conjunto de prueba, en general, el modelo con mayor área bajo la curva presentará los mejores resultados de clasificación.
linear <- svm(Purchase ~ ., data = oj_train, kernel = "linear", probability = TRUE, cost = tune_oj$best.model$cost)
pred_linear <- predict(linear, newdata = oj_test, probability = TRUE)
radial <- svm(Purchase ~ ., data = oj_train, kernel = "radial", probability = TRUE, cost = tune_oj_radial$best.model$cost)
pred_radial <- predict(radial, newdata = oj_test, probability = TRUE)
poly <- svm(Purchase ~ ., data = oj_train, kernel = "polynomial", probability = TRUE, cost = tune_oj_poly$best.model$cost)
pred_poly <- predict(poly, newdata = oj_test, probability = TRUE)
roc_linear <- tibble(
estimate = pred_linear,
truth = oj_test$Purchase,
prob = attr(pred_linear,"probabilities")[, 2],
Modelo = "Lineal")
roc_radial <- tibble(
estimate = pred_radial,
truth = oj_test$Purchase,
prob = attr(pred_radial,"probabilities")[, 2],
Modelo = "Radial")
roc_poly <- tibble(
estimate = pred_poly,
truth = oj_test$Purchase,
prob = attr(pred_poly,"probabilities")[, 2],
Modelo = "Polinomial")Con la utilización del argumento probability = TRUE se le indica al método svm que calcule las probabilidades necesarias para la curva ROC.
roc_merge <- rbind(roc_linear, roc_radial, roc_poly)
roc_merge %>%
group_by(Modelo) %>%
roc_curve(truth, prob) %>%
autoplot() +
labs(title = "Curvas ROC de los modelos SVM con diferente Kernel")En general se identifica que la curva del modelo Polinomial se enceuntra por debajo de las demás. Se puede notar que los modelos entrenados tienen una buena calidad pues las curvas ROC estan alejadas de la recta dada por la función \(f(x) = x\).
Con el objetivo de determinar una diferencia entre los modelos, veamos el AUC para cada uno de ellos.
## # A tibble: 3 x 4
## Modelo .metric .estimator .estimate
## <chr> <chr> <chr> <dbl>
## 1 Lineal roc_auc binary 0.906
## 2 Polinomial roc_auc binary 0.888
## 3 Radial roc_auc binary 0.882
El AUC del modelo lineal es ligeramente mayor que el de los demás modelos, lo cual indica que la suposición de que la complejidad/flexibilidad añadida por los kernel radial y polinomial al haber gran cantidad de variables busca patrones en los datos de entrenamiento que posiblemente no son generales para para el conjunto de datos completo, generando un ligero sobreajuste.
Conclusiones
Durante el desarrollo de este trabajo se aplicaron diferentes técnicas de aprendizaje estadístico (Regresión logística, Regresión lineal, Regresión lineal múltiple, Stepwise regression, Árboles de decisión, bosques aleatorios, Boosting, Bagging, Máquinas de soporte vectorial) sobre conjuntos de datos de la vida real y también generados sintéticamente; con el objetivo de entender el principio teórico de los algoritmos y además afrontar las dificultades técnicas que puede llevar su implementación.
Dada la extensión de este trabajo (Casi 3000 líneas de código y comentarios), aprendimos a trabajar en equipo de manera muy integrada utilizando control de versiones mediante Git y verificando el trabajo realizado, afrontamos además las dificultades de compilar un informe de gran extensión y los errores que pudieron surgir en el camino.
En este trabajo hubo ejercicios en los que se realizó análisis descriptivo exhaustivo, gracias a esto identificamos que éste es un paso esencial que puede dar muchas entradas para la selección de un modelo correcto y las variables que se deberían incorporar o no. Además, aprendimos a poner los diferentes modelos a prueba y a identificar cuál es el más adecuado utilizando medidas de error y visualizaciones correctas.
En muchos de estos ejercicios identificamos que a veces los métodos más clásicos como la regresión lineal pueden entregar resultados mejores que algoritmos de alta complejidad, por lo que siempre será de nuestro interés verificar las circunstancias bajo las que se decide utilizar una técnica más compleja y evaluar que el costo añadido de cómputo, abstracción y programación valga la pena.